Repository: bumps/bumps (language: python)

Issue:
restructure parameter and expression handling

minimal example of model constraints

Check that fit weighting works
Looking at the code I don't think weighted fitting will work:
https://github.com/bumps/bumps/blob/a1d811a7d1ffd3f07800659edb0bbbd974ab48c1/bumps/fitproblem.py#L683
It should be:
```
return sum(w*f.model_nllf() for w, f in zip(self.weights, self.models))
```

There may be additional issues. Test cases are needed to maintain this feature.
covariance calculation is incorrect for 2D data
After fitting `doc/examples/peaks/model.py` with amoeba or DE a Jacobian is computed at the fitted optimum. To build the Jacobian a forward difference calculation is performed, with the resulting residuals stacked as rows in a matrix then transposed. Because the residuals are 2D the resulting matrix has the wrong shape, and the covariance has the wrong number of columns.

The fix is to flatten the residuals array before stacking the matrix:
https://github.com/bumps/bumps/blob/137cd261755b47824258a43dc11c1895ff5c418b/bumps/lsqerror.py#L77-L79

With `f` returning a row vector this will fix the forward difference calculation here:
https://github.com/bumps/bumps/blob/137cd261755b47824258a43dc11c1895ff5c418b/bumps/lsqerror.py#L100

as well as the central difference calculation here:
https://github.com/bumps/bumps/blob/137cd261755b47824258a43dc11c1895ff5c418b/bumps/lsqerror.py#L117

The Hessian is fine since it uses the scalar function nllf. 



Discussion / hints:

Fix patch:
diff --git a/bumps/bounds.py b/bumps/bounds.py
index bf6f92d2..d7b6cdfa 100644
--- a/bumps/bounds.py
+++ b/bumps/bounds.py
@@ -66,8 +66,24 @@
     # failure if it doesn't exist.
     pass
 
+from .util import field, schema, Optional, Any, Union, Dict, Callable, Literal, Tuple, List, Literal
 
-def pm(v, plus, minus=None, limits=None):
+LimitValue = Union[float, Literal["-inf", "inf"]]
+LimitsType = Tuple[Union[float, Literal["-inf"]], Union[float, Literal["inf"]]]
+
+# TODO: should we use this in the bounds limits?
+# @schema()
+# class ExtendedFloat(float):
+#     __root__: Union[float, Literal["inf", "-inf"]]
+#     def __new__(cls, *args, **kw):
+#         return super().__new__(cls, *args)
+#     def __init__(self, __root__=None):
+#         pass
+#     def __repr__(self):
+#         return float.__repr__(self)
+
+
+def pm(v, plus, minus=None, limits: Optional[LimitsType]=None):
     """
     Return the tuple (~v-dv,~v+dv), where ~expr is a 'nice' number near to
     to the value of expr.  For example::
@@ -151,7 +167,7 @@ def nice_range(bounds):
         return bounds
 
 
-def init_bounds(v):
+def init_bounds(v) -> 'Bounds':
     """
     Returns a bounds object of the appropriate type given the arguments.
 
@@ -186,7 +202,7 @@ def init_bounds(v):
         return Bounded(lo, hi)
 
 
-class Bounds(object):
+class Bounds:
     """
     Bounds abstract base class.
 
@@ -201,9 +217,12 @@ class Bounds(object):
     is being optimized is also a probability, then this is an easy way to
     incorporate information from other sorts of measurements into the model.
     """
-    limits = (-inf, inf)
     # TODO: need derivatives wrt bounds transforms
 
+    @property
+    def limits(self):
+        return (-inf, inf)
+
     def get01(self, x):
         """
         Convert value into [0,1] for optimizers which are bounds constrained.
@@ -281,6 +300,19 @@ def __str__(self):
         limits = tuple(num_format(v) for v in self.limits)
         return "(%s,%s)" % limits
 
+    def satisfied(self, v) -> bool:
+        lo, hi = self.limits
+        return v >= lo and v <= hi
+
+    def penalty(self, v) -> float:
+        """
+        return a (differentiable) nonzero value when outside the bounds
+        """
+        lo, hi = self.limits
+        dlo = 0. if v >= lo else abs(v - lo)
+        dhi = 0. if v <= hi else abs(v - hi)
+        return dlo + dhi
+
     def to_dict(self):
         return dict(
             type=type(self).__name__,
@@ -302,6 +334,7 @@ def num_format(v):
         return "NaN"
 
 
+@schema()
 class Unbounded(Bounds):
     """
     Unbounded parameter.
@@ -312,6 +345,10 @@ class Unbounded(Bounds):
     log likelihood of P is inf everywhere.  A value inf will interfere
     with optimization routines, and so we instead choose P == 1 everywhere.
     """
+    type = "Unbounded"
+
+    def __init__(self, *args, **kw):
+        pass
 
     def random(self, n=1, target=1.0):
         scale = target + (target == 0.)
@@ -336,6 +373,7 @@ def putfull(self, v):
         return v
 
 
+@schema()
 class BoundedBelow(Bounds):
     """
     Semidefinite range bounded below.
@@ -354,27 +392,32 @@ class BoundedBelow(Bounds):
     is indistinguishable from values outside the range.  Instead we say
     that P = 1 in range, and 0 outside.
     """
+    base: float
+    type = "BoundedBelow"
 
     def __init__(self, base):
-        self.limits = (base, inf)
-        self._base = base
+        self.base = base
+
+    @property
+    def limits(self):
+        return (self.base, inf)
 
     def start_value(self):
-        return self._base + 1
+        return self.base + 1
 
-    def random(self, n=1, target=1.):
-        target = max(abs(target), abs(self._base))
-        scale = target + (target == 0.)
-        return self._base + abs(RNG.randn(n)*scale)
+    def random(self, n=1, target:float=1.):
+        target = max(abs(target), abs(self.base))
+        scale = target + float(target == 0.)
+        return self.base + abs(RNG.randn(n)*scale)
 
     def nllf(self, value):
-        return 0 if value >= self._base else inf
+        return 0 if value >= self.base else inf
 
     def residual(self, value):
-        return 0 if value >= self._base else -4
+        return 0 if value >= self.base else -4
 
     def get01(self, x):
-        m, e = math.frexp(x - self._base)
+        m, e = math.frexp(x - self.base)
         if m >= 0 and e <= _E_MAX:
             v = (e + m) / (2. * _E_MAX)
             return v
@@ -385,18 +428,19 @@ def put01(self, v):
         v = v * 2 * _E_MAX
         e = int(v)
         m = v - e
-        x = math.ldexp(m, e) + self._base
+        x = math.ldexp(m, e) + self.base
         return x
 
     def getfull(self, x):
-        v = x - self._base
+        v = x - self.base
         return v if v >= 1 else 2 - 1. / v
 
     def putfull(self, v):
         x = v if v >= 1 else 1. / (2 - v)
-        return x + self._base
+        return x + self.base
 
 
+@schema()
 class BoundedAbove(Bounds):
     """
     Semidefinite range bounded above.
@@ -413,27 +457,31 @@ class BoundedAbove(Bounds):
     is indistinguishable from values outside the range.  Instead we say
     that P = 1 in range, and 0 outside.
     """
+    base: float
 
     def __init__(self, base):
-        self.limits = (-inf, base)
-        self._base = base
+        self.base = base
+
+    @property
+    def limits(self):
+        return (-inf, self.base)
 
     def start_value(self):
-        return self._base - 1
+        return self.base - 1
 
-    def random(self, n=1, target=1.0):
-        target = max(abs(self._base), abs(target))
-        scale = target + (target == 0.)
-        return self._base - abs(RNG.randn(n)*scale)
+    def random(self, n=1, target:float=1.0):
+        target = max(abs(self.base), abs(target))
+        scale = target + float(target == 0.)
+        return self.base - abs(RNG.randn(n)*scale)
 
     def nllf(self, value):
-        return 0 if value <= self._base else inf
+        return 0 if value <= self.base else inf
 
     def residual(self, value):
-        return 0 if value <= self._base else 4
+        return 0 if value <= self.base else 4
 
     def get01(self, x):
-        m, e = math.frexp(self._base - x)
+        m, e = math.frexp(self.base - x)
         if m >= 0 and e <= _E_MAX:
             v = (e + m) / (2. * _E_MAX)
             return 1 - v
@@ -444,18 +492,19 @@ def put01(self, v):
         v = (1 - v) * 2 * _E_MAX
         e = int(v)
         m = v - e
-        x = -(math.ldexp(m, e) - self._base)
+        x = -(math.ldexp(m, e) - self.base)
         return x
 
     def getfull(self, x):
-        v = x - self._base
+        v = x - self.base
         return v if v <= -1 else -2 - 1. / v
 
     def putfull(self, v):
         x = v if v <= -1 else -1. / (v + 2)
-        return x + self._base
+        return x + self.base
 
 
+@schema()
 class Bounded(Bounds):
     """
     Bounded range.
@@ -468,24 +517,33 @@ class Bounded(Bounds):
     and for a more natural mapping between nllf and chisq, we instead
     set the probability to 0.  This choice will not affect the fits.
     """
+    lo: float = field(metadata={"description": "lower end of bounds"})
+    hi: float = field(metadata={"description": "upper end of bounds"})
+
+    # @classmethod
+    # def from_dict(cls, limits=None):
+    #     lo, hi = limits
+    #     return cls(lo, hi)
 
     def __init__(self, lo, hi):
-        self.limits = (lo, hi)
+        self.lo = lo
+        self.hi = hi
         self._nllf_scale = log(hi - lo)
 
+    @property
+    def limits(self):
+        return (self.lo, self.hi)
+
     def random(self, n=1, target=1.0):
-        lo, hi = self.limits
         #print("= uniform",lo,hi)
-        return RNG.uniform(lo, hi, size=n)
+        return RNG.uniform(self.lo, self.hi, size=n)
 
     def nllf(self, value):
-        lo, hi = self.limits
-        return 0 if lo <= value <= hi else inf
+        return 0 if self.lo <= value <= self.hi else inf
         # return self._nllf_scale if lo<=value<=hi else inf
 
     def residual(self, value):
-        lo, hi = self.limits
-        return -4 if lo > value else (4 if hi < value else 0)
+        return -4 if self.lo > value else (4 if self.hi < value else 0)
 
     def get01(self, x):
         lo, hi = self.limits
@@ -502,6 +560,7 @@ def putfull(self, v):
         return self.put01(_get01_inf(v))
 
 
+
 class Distribution(Bounds):
     """
     Parameter is pulled from a distribution.
@@ -510,6 +569,8 @@ class Distribution(Bounds):
     In particular, it should define methods rvs, nnlf, cdf and ppf and
     attributes args and dist.name.
     """
+    #dist: DistProtocol = None
+    #dist: Any = None
 
     def __init__(self, dist):
         self.dist = dist
@@ -556,6 +617,7 @@ def to_dict(self):
             )
 
 
+@schema()
 class Normal(Distribution):
     """
     Parameter is pulled from a normal distribution.
@@ -568,10 +630,12 @@ class Normal(Distribution):
     *mean* is the expected value of the parameter and *std* is the 1-sigma
     standard deviation.
     """
+    mean: float = 0.0
+    std: float = 1.0
 
-    def __init__(self, mean=0, std=1):
+    def __init__(self, mean:float=0, std:float=1):
         Distribution.__init__(self, normal_distribution(mean, std))
-        self._nllf_scale = log(sqrt(2 * pi * std ** 2))
+        self._nllf_scale = log(2 * pi * std ** 2)/2
 
     def nllf(self, value):
         # P(v) = exp(-0.5*(v-mean)**2/std**2)/sqrt(2*pi*std**2)
@@ -592,18 +656,34 @@ def __setstate__(self, state):
         self.__init__(mean=mean, std=std)
 
 
+@schema()
 class BoundedNormal(Bounds):
     """
     truncated normal bounds
     """
-
-    def __init__(self, mean=0, std=1, limits=(-inf, inf)):
-        self.limits = limits
+    mean: float = 0.0
+    std: float = 1.0
+    lo: Union[float, Literal["-inf"]]
+    hi: Union[float, Literal["inf"]]
+
+    def __init__(self, mean:float=0, std:float=1, limits=(-inf, inf), hi="inf", lo="-inf"):
+        if limits is not None:
+            # for backward compatibility:
+            lo, hi = limits
+        limits = (
+            -inf if lo is None else float(lo),
+            inf if hi is None else float(hi)
+        )
+        self.lo, self.hi = limits
         self.mean, self.std = mean, std
 
         self._left = normal_distribution.cdf((limits[0]-mean)/std)
         self._delta = normal_distribution.cdf((limits[1]-mean)/std) - self._left
-        self._nllf_scale = log(sqrt(2 * pi * std ** 2)) + log(self._delta)
+        self._nllf_scale = log(2 * pi * std ** 2)/2 + log(self._delta)
+
+    @property
+    def limits(self):
+        return (self.lo, self.hi)
 
     def get01(self, x):
         """
@@ -680,6 +760,7 @@ def __str__(self):
         return "(%s,%s), norm(%s,%s)" % tuple(num_format(v) for v in vals)
 
 
+@schema()
 class SoftBounded(Bounds):
     """
     Parameter is pulled from a stretched normal distribution.
@@ -695,41 +776,48 @@ class SoftBounded(Bounds):
     into the range [0,1] for each parameter we don't need to use soft
     constraints, and this acts just like the rectangular distribution.
     """
+    lo: float = 0.0
+    hi: float = 1.0
+    std: float = 1.0
 
-    def __init__(self, lo, hi, std=None):
-        self._lo, self._hi, self._std = lo, hi, std
+    def __init__(self, lo, hi, std=1.0):
+        self.lo, self.hi, self.std = lo, hi, std
         self._nllf_scale = log(hi - lo + sqrt(2 * pi * std))
 
+    @property
+    def limits(self):
+        return (self.lo, self.hi)
+
     def random(self, n=1, target=1.0):
-        return RNG.uniform(self._lo, self._hi, size=n)
+        return RNG.uniform(self.lo, self.hi, size=n)
 
     def nllf(self, value):
         # To turn f(x) = 1 if x in [lo,hi] else G(tail)
         # into a probability p, we need to normalize by \int{f(x)dx},
         # which is just hi-lo + sqrt(2*pi*std**2).
-        if value < self._lo:
-            z = self._lo - value
-        elif value > self._hi:
-            z = value - self._hi
+        if value < self.lo:
+            z = self.lo - value
+        elif value > self.hi:
+            z = value - self.hi
         else:
             z = 0
-        return (z / self._std) ** 2 / 2 + self._nllf_scale
+        return (z / self.std) ** 2 / 2 + self._nllf_scale
 
     def residual(self, value):
-        if value < self._lo:
-            z = self._lo - value
-        elif value > self._hi:
-            z = value - self._hi
+        if value < self.lo:
+            z = self.lo - value
+        elif value > self.hi:
+            z = value - self.hi
         else:
             z = 0
-        return z / self._std
+        return z / self.std
 
     def get01(self, x):
-        v = float(x - self._lo) / (self._hi - self._lo)
+        v = float(x - self.lo) / (self.hi - self.lo)
         return v if 0 <= v <= 1 else (0 if v < 0 else 1)
 
     def put01(self, v):
-        return v * (self._hi - self._lo) + self._lo
+        return v * (self.hi - self.lo) + self.lo
 
     def getfull(self, x):
         return x
@@ -738,7 +826,7 @@ def putfull(self, v):
         return v
 
     def __str__(self):
-        return "box_norm(%g,%g,sigma=%g)" % (self._lo, self._hi, self._std)
+        return "box_norm(%g,%g,sigma=%g)" % (self.lo, self.hi, self.std)
 
 
 _E_MIN = -1023
@@ -784,3 +872,5 @@ def _put01_inf(v):
     x = math.ldexp(s * m, e + _E_MIN)
     # print "< x,e,m,s,v",x,e+_e_min,s*m,s,v
     return x
+
+BoundsType = Union[Unbounded, Bounded, BoundedAbove, BoundedBelow, BoundedNormal, SoftBounded, Normal]
\ No newline at end of file
diff --git a/bumps/dream/corrplot.py b/bumps/dream/corrplot.py
index ccd42cf4..439f4e21 100644
--- a/bumps/dream/corrplot.py
+++ b/bumps/dream/corrplot.py
@@ -49,14 +49,15 @@ def __getitem__(self, i, j):
         """
         return self.hists[i, j]
 
-    def plot(self, title=None):
+    def plot(self, title=None, fig=None):
         """
         Plot the correlation histograms on the specified figure
         """
         import pylab
 
-        pylab.clf()
-        fig = pylab.gcf()
+        if fig is None:
+            pylab.clf()
+            fig = pylab.gcf()
         if title is not None:
             fig.text(0.5, 0.95, title,
                      horizontalalignment='center',
diff --git a/bumps/dream/entropy.py b/bumps/dream/entropy.py
index c5b52136..f965f781 100644
--- a/bumps/dream/entropy.py
+++ b/bumps/dream/entropy.py
@@ -859,8 +859,14 @@ def rvs(self, size=1):
         return self.dist.rvs(size=size)
 
     def entropy(self, N=10000):
+        # CRUFT scipy==1.10.0: scipy.stats briefly removed the dist.cov attribute.
+        if hasattr(self.dist, 'cov'):
+            cov = self.dist.cov
+        else:
+            cov = self.dist.cov_object.covariance
+
         with np.errstate(divide='ignore'):
-            return 0.5*log(np.linalg.det((2*pi*np.e)*self.dist.cov))
+            return 0.5*log(np.linalg.det((2*pi*np.e)*cov))
 
 class GaussianMixture:
     def __init__(self, w, mu=None, sigma=None):
diff --git a/bumps/dream/varplot.py b/bumps/dream/varplot.py
index d27038dc..17941dc4 100644
--- a/bumps/dream/varplot.py
+++ b/bumps/dream/varplot.py
@@ -34,13 +34,14 @@ def var_plot_size(n):
     figheight = (TILE_H+V_SPACE)*nrow + T_MARGIN + B_MARGIN
     return figwidth, figheight
 
-def _make_var_axes(n):
+def _make_var_axes(n, fig=None):
     """
     Build a figure with one axis per parameter,
     and one axis (the last one) to contain the colorbar.
     Use to make the vars histogram figure.
     """
-    fig = plt.gcf()
+    if fig is None:
+        fig = plt.gcf()
     fig.clf()
     total_width, total_height = fig.get_size_inches()
 
@@ -96,36 +97,39 @@ def tile_axes_square(n):
     return cols, rows
 
 
-def plot_vars(draw, all_vstats, **kw):
+def plot_vars(draw, all_vstats, fig=None, **kw):
     n = len(all_vstats)
-    fig = _make_var_axes(n)
-    cbar = _make_fig_colorbar(draw.logp)
+    fig = _make_var_axes(n, fig=fig)
+    cbar = _make_fig_colorbar(draw.logp, fig=fig)
     for k, vstats in enumerate(all_vstats):
-        fig.sca(fig.axes[k])
-        plot_var(draw, vstats, k, cbar, **kw)
+        axes = fig.axes[k]
+        plot_var(draw, vstats, k, cbar, axes=axes, **kw)
+        fig.canvas.draw()
 
 
-def plot_var(draw, vstats, var, cbar, nbins=30):
+def plot_var(draw, vstats, var, cbar, nbins=30, axes=None):
     values = draw.points[:, var].flatten()
     bin_range = vstats.p95_range
     #bin_range = np.min(values), np.max(values)
+    import pylab
+    if axes is None:
+        axes = pylab.gca()
+
     _make_logp_histogram(values, draw.logp, nbins, bin_range,
-                         draw.weights, cbar)
-    _decorate_histogram(vstats)
+                         draw.weights, cbar, axes)
+    _decorate_histogram(vstats, axes)
 
 
-def _decorate_histogram(vstats):
-    import pylab
+def _decorate_histogram(vstats, axes):
     from matplotlib.transforms import blended_transform_factory as blend
 
     l95, h95 = vstats.p95_range
     l68, h68 = vstats.p68_range
 
     # Shade things inside 1-sigma
-    pylab.axvspan(l68, h68, color='gold', alpha=0.5, zorder=-1)
+    axes.axvspan(l68, h68, color='gold', alpha=0.5, zorder=-1)
     # build transform with x=data, y=axes(0,1)
-    ax = pylab.gca()
-    transform = blend(ax.transData, ax.transAxes)
+    transform = blend(axes.transData, axes.transAxes)
 
     def marker(symbol, position):
         if position < l95:
@@ -134,23 +138,23 @@ def marker(symbol, position):
             symbol, position, ha = '>'+symbol, h95, 'right'
         else:
             symbol, position, ha = symbol, position, 'center'
-        pylab.text(position, 0.95, symbol, va='top', ha=ha,
+        axes.text(position, 0.95, symbol, va='top', ha=ha,
                    transform=transform, zorder=3, color='g')
-        #pylab.axvline(v)
+        #axes.axvline(v)
 
     marker('|', vstats.median)
     marker('E', vstats.mean)
     marker('*', vstats.best)
 
-    pylab.text(0.01, 0.95, vstats.label, zorder=2,
+    axes.text(0.01, 0.95, vstats.label, zorder=2,
                backgroundcolor=(1, 1, 0, 0.2),
                verticalalignment='top',
                horizontalalignment='left',
-               transform=pylab.gca().transAxes)
-    ax.set_yticklabels([])
+               transform=axes.transAxes)
+    axes.set_yticklabels([])
 
 
-def _make_fig_colorbar(logp):
+def _make_fig_colorbar(logp, fig=None):
     import matplotlib as mpl
     import pylab
 
@@ -164,7 +168,8 @@ def _make_fig_colorbar(logp):
     # Option 3: full range
     #vmin,vmax = -max(logp),-min(logp)
 
-    fig = pylab.gcf()
+    if fig is None:
+        fig = pylab.gcf()
     ax = fig.axes[-1]
     cmap = mpl.cm.copper
 
@@ -182,7 +187,10 @@ def __init__(self, low, high):
             self.delta = high - low
 
         def __call__(self, x, pos=None):
-            return format_value(x, self.delta)
+            # TODO: where did format_value come from?
+            # it does not exist anywhere in the project.
+            # return format_value(x, self.delta)
+            return '{:.3G}'.format(x)
 
     ticks = ()  #(vmin, vmax)
     formatter = MinDigitsFormatter(vmin, vmax)
@@ -203,7 +211,7 @@ def __call__(self, x, pos=None):
     return cbar
 
 
-def _make_logp_histogram(values, logp, nbins, ci, weights, cbar):
+def _make_logp_histogram(values, logp, nbins, ci, weights, cbar, axes):
     from numpy import (ones_like, searchsorted, linspace, cumsum, diff,
                        unique, argsort, array, hstack, exp)
     if weights is None:
@@ -218,7 +226,6 @@ def _make_logp_histogram(values, logp, nbins, ci, weights, cbar):
     #weightsum = cumsum(weights)
     #heights = diff(weightsum[idx])/weightsum[-1]  # normalized weights
 
-    import pylab
     edgecolors = None
     cmap = cbar.cmap
     cmap_edges = linspace(0, 1, cmap.N+1)[1:-1]
@@ -274,7 +281,7 @@ def _make_logp_histogram(values, logp, nbins, ci, weights, cbar):
         tops = unique(hstack((change_point, len(pv)-1)))
         y = hstack((0, y_top[tops]))
         z = pv[tops][:, None]
-        pylab.pcolormesh(
+        axes.pcolormesh(
             x, y, z, norm=cbar.norm, cmap=cmap, edgecolors=edgecolors)
 
         # centerpoint, histogram height, maximum likelihood for each bin
@@ -291,7 +298,7 @@ def _make_logp_histogram(values, logp, nbins, ci, weights, cbar):
     ml_peak = np.max(maxlikelihood)
     if ml_peak > hist_peak*1.3:
         maxlikelihood *= hist_peak*1.3/ml_peak
-    pylab.plot(centers, maxlikelihood, '-g')
+    axes.plot(centers, maxlikelihood, '-g')
 
     ## plot marginal gaussian approximation along with histogram
     #def G(x, mean, std):
diff --git a/bumps/dream/views.py b/bumps/dream/views.py
index 7aded468..a6f1456a 100644
--- a/bumps/dream/views.py
+++ b/bumps/dream/views.py
@@ -84,9 +84,9 @@ def plot_all(state, portion=1.0, figfile=None):
         if figfile is not None:
             savefig(figfile+"-parcor"+figext)
 
-def plot_corrmatrix(draw, nbins=50):
+def plot_corrmatrix(draw, nbins=50, fig=None):
     c = corrplot.Corr2d(draw.points.T, bins=nbins, labels=draw.labels)
-    c.plot()
+    c.plot(fig=fig)
     #print "Correlation matrix\n",c.R()
 
 
diff --git a/bumps/errplot.py b/bumps/errplot.py
index 485eb7c7..80ca1438 100644
--- a/bumps/errplot.py
+++ b/bumps/errplot.py
@@ -123,10 +123,10 @@ def calc_errors(problem, points):
     return ret
 
 
-def show_errors(errs):
+def show_errors(errs, fig=None):
     """
     Display the confidence regions returned by :func:`calc_errors`.
 
     The content of *errs* depends on the active plugin.
     """
-    return plugin.show_errors(errs)
+    return plugin.show_errors(errs, fig=fig)
diff --git a/bumps/fitproblem.py b/bumps/fitproblem.py
index 5e430ba2..ca9a9679 100644
--- a/bumps/fitproblem.py
+++ b/bumps/fitproblem.py
@@ -49,31 +49,41 @@
 # to avoid this issue.
 from __future__ import division, with_statement
 
-__all__ = ['Fitness', 'FitProblem', 'load_problem',
-           'BaseFitProblem', 'MultiFitProblem']
+__all__ = ['Fitness', 'FitProblem', 'load_problem']
 
 import sys
 import os
 import traceback
 import logging
 import warnings
+import builtins
+
+from typing import Optional
 
 import numpy as np
-from numpy import inf, isnan, NaN
+from numpy import inf, isnan, NaN, isfinite
 
 from . import parameter, bounds as mbounds
-from .parameter import to_dict
+from .parameter import to_dict, Parameter, Expression, Variable
 from .formatnum import format_uncertainty
 from . import util
 
-# Abstract base class
-class Fitness(object):
+if 'SPHINXOPTS' in os.environ:
+    # export this only when building docs
+    __all__.append('FitProblemSchema')
+
+# Abstract base class:
+# can use "isinstance" to check if a class implements the protocol
+@util.runtime_checkable
+@util.schema(init=False)
+class Fitness(util.Protocol):
     """
     Manage parameters, data, and theory function evaluation.
 
     See :ref:`fitness` for a detailed explanation.
     """
-    def parameters(self):
+
+    def parameters(self) -> util.List[Parameter]:
         """
         return the parameters in the model.
 
@@ -143,22 +153,42 @@ def plot(self, view='linear'):
         pass
 
 
-def no_constraints():
+def no_constraints() -> float:
     """default constraints function for FitProblem"""
     return 0
 
+def chisq_str(fitness: Fitness) -> str:
+    """
+    Return a string representing the chisq equivalent of the nllf.
+    If the model has strictly gaussian independent uncertainties then the
+    negative log likelihood function will return 0.5*sum(residuals**2),
+    which is 1/2*chisq.  Since we are printing normalized chisq, we
+    multiply the model nllf by 2/DOF before displaying the value.  This
+    is different from the problem nllf function, which includes the
+    cost of the prior parameters and the cost of the penalty constraints
+    in the total nllf.  The constraint value is displayed separately.
+    """
 
-# TODO: refactor FitProblem definition
-# deprecate the direct use of MultiFitProblem
-def FitProblem(*args, **kw):
+    chisq = fitness.nllf() / fitness.numpoints()
+    text = format_uncertainty(chisq, None)
+    return text
+
+def show_parameters(fitness: Fitness, subs: Optional[util.Dict[util.Any, Parameter]]=None):
+    """Print the available parameters to the console as a tree."""
+    print(parameter.format(fitness.parameters(), freevars=subs))
+    print("[chisq=%s, nllf=%g]" % (chisq_str(fitness), fitness.nllf()))
+
+@util.schema(classname="FitProblem", init=False, eq=False)
+class FitProblemSchema:
+    name: util.Optional[str] = None
+    models: util.List[Fitness]
+    freevars: util.Optional[parameter.FreeVariables] = None
+    weights: util.Union[util.List[float], util.Literal[None]] = None
+    constraints: util.Optional[util.Sequence[parameter.Constraint]] = None
+    penalty_nllf: util.Union[float, util.Literal["inf"]] = "inf"
+
+class FitProblem(FitProblemSchema):
     r"""
-    Return a fit problem instance for the fitness function(s).
-
-    For an individual model:
-
-        *fitness* is a :class:`Fitness` instance.
-
-    For a set of models:
 
         *models* is a sequence of :class:`Fitness` instances.
 
@@ -177,7 +207,10 @@ def FitProblem(*args, **kw):
 
         *name* name of the problem
 
-        *constraints* is a function which returns the negative log likelihood
+        *constraints* is a list of Constraint objects, which have a method to
+        calculate the nllf for that constraint. 
+        Also supports an alternate form which cannot be serialized:
+        A function which returns the negative log likelihood
         of seeing the parameters independent from the fitness function.  Use
         this for example to check for feasible regions of the search space, or
         to add constraints that cannot be easily calculated per parameter.
@@ -185,11 +218,8 @@ def FitProblem(*args, **kw):
         the feasible region so that the fit will be directed toward feasible
         values.
 
-        *soft_limit* is the constraints function cutoff, beyond which the
-        *penalty_nllf* will be used and *fitness* nllf will not be calculated.
-
         *penalty_nllf* is the nllf to use for *fitness* when *constraints*
-        is greater than *soft_limit*.
+        or model parameter bounds are not satisfied.
 
     Total nllf is the sum of the parameter nllf, the constraints nllf and the
     depending on whether constraints is greater than soft_limit, either the
@@ -197,119 +227,100 @@ def FitProblem(*args, **kw):
 
     New in 0.9.0: weights are now squared when computing the sum rather than
     linear.
-    """
-    if len(args) > 0:
-        if isinstance(args[0], (list, tuple)):
-            return MultiFitProblem(args[0], *args[1:], **kw)
-        else:
-            return BaseFitProblem(*args, **kw)
-    else:
-        if 'fitness' in kw:
-            return BaseFitProblem(*args, **kw)
-        else:
-            return MultiFitProblem(*args, **kw)
+    """  
 
+    _constraints_function: util.Callable[..., float]
+    _parameters: util.List[Parameter]
+    _parameters_by_id: util.Dict[str, Parameter]
+    freevars: parameter.FreeVariables
+    # _all_constraints: util.List[util.Union[Parameter, Expression]]
 
-class BaseFitProblem(object):
-    """
-    See :func:`FitProblem`
-    """
-    def __init__(self, fitness, name=None, constraints=None,
-                 penalty_nllf=np.inf, soft_limit=np.inf, partial=False):
-        self.constraints = constraints
-        self.fitness = fitness
-        self.partial = partial
-        if name is not None:
-            self.name = name
+    def __init__(self, models: util.Union[Fitness, util.List[Fitness]], weights=None, name=None,
+                 constraints=None,
+                 penalty_nllf="inf",
+                 freevars=None,
+                 soft_limit:Optional[float]=None,  # TODO: deprecate
+                 ):
+        if not isinstance(models, (list, tuple)):
+            models = [models]
+        if callable(constraints):
+            warnings.warn("convert constraints functions to constraint expressions", DeprecationWarning)
+            self._constraints_function = constraints
+            self.constraints = []
         else:
-            try:
-                self.name = fitness.name
-            except AttributeError:
-                self.name = 'FitProblem'
+            self._constraints_function = self._null_constraints_function
+            # TODO: do we want to allow "constraints=a<b" or do we require a sequence "constraints=[a<b]"?
+            self.constraints = constraints if constraints is not None else []
+        if freevars is None:
+            names = ["M%d" % i for i, _ in enumerate(models)]
+            freevars = parameter.FreeVariables(names=names)
+        self.freevars = freevars
+        self._models = models
+        if weights is None:
+            weights = [1.0 for _ in models]
+        self.weights = weights
+        self.penalty_nllf = float(penalty_nllf)
+        self.set_active_model(0)  # Set the active model to model 0
+        self.model_reset() # sets self._all_constraints
+        self.name = name
 
-        self.soft_limit = soft_limit
-        self.penalty_nllf = penalty_nllf
-        self.model_reset()
+    @staticmethod
+    def _null_constraints_function():
+        return 0.
 
-    # noinspection PyAttributeOutsideInit
-    def model_reset(self):
-        """
-        Prepare for the fit.
+    @property
+    def fitness(self):
+        warnings.warn('Deprecated: use of problem.fitness will be removed at some point')
+        if len(self._models) == 1:
+            return self._models[0]
+        raise ValueError('problem.fitness is not defined')
 
-        This sets the parameters and the bounds properties that the
-        solver is expecting from the fittable object.  We also compute
-        the degrees of freedom so that we can return a normalized fit
-        likelihood.
+    # TODO: make this @property\ndef models(self): ...
+    @property
+    def models(self):
+        """Iterate over models, with free parameters set from model values"""
+        try:
+            for i, f in enumerate(self._models):
+                self.freevars.set_model(i)
+                yield f
+        finally:
+            # Restore the active model after cycling, even if interrupted
+            self.freevars.set_model(self._active_model_index)
 
-        If the set of fit parameters changes, then model_reset must
-        be called.
-        """
-        # print self.model_parameters()
-        all_parameters = parameter.unique(self.model_parameters())
-        # print "all_parameters",all_parameters
-        self._parameters = parameter.varying(all_parameters)
-        # print "varying",self._parameters
-        self.bounded = [p for p in all_parameters
-                        if not isinstance(p.bounds, mbounds.Unbounded)]
-        self.dof = self.model_points()
-        if not self.partial:
-            self.dof -= len(self._parameters)
-        if self.dof <= 0:
-            raise ValueError("Need more data points than fitting parameters")
-        #self.constraints = pars.constraints()
+    # noinspection PyAttributeOutsideInit
+    def set_active_model(self, i):
+        """Use free parameters from model *i*"""
+        self._active_model_index = i
+        self.active_model = self._models[i]
+        self.freevars.set_model(i)
 
     def model_parameters(self):
-        """
-        Parameters associated with the model.
-        """
-        return self.fitness.parameters()
+        """Return parameters from all models"""
+        pars = {}
+        pars['models'] = [f.parameters() for f in self.models]
+        free = self.freevars.parameters()
+        if free:
+            pars['freevars'] = free
+        return pars
 
     def to_dict(self):
         return {
             'type': type(self).__name__,
             'name': self.name,
-            'fitness': to_dict(self.fitness),
-            'partial': self.partial,
-            'soft_limit': self.soft_limit,
+            'models': to_dict(self._models),
+            'weights': self.weights,
             'penalty_nllf': self.penalty_nllf,
             # TODO: constraints may be a function.
             'constraints': to_dict(self.constraints),
+            'freevars': to_dict(self.freevars),
         }
 
-    def model_points(self):
-        """
-        Number of data points associated with the model.
-        """
-        return self.fitness.numpoints()
-
-    def model_update(self):
-        """
-        Update the model according to the changed parameters.
-        """
-        if hasattr(self.fitness, 'update'):
-            self.fitness.update()
-
-    def model_nllf(self):
-        """
-        Negative log likelihood of seeing data given model.
-        """
-        return self.fitness.nllf()
-
-    def simulate_data(self, noise=None):
-        """Simulate data with added noise"""
-        self.fitness.simulate_data(noise=noise)
-
-    def resynth_data(self):
-        """Resynthesize data with noise from the uncertainty estimates."""
-        self.fitness.resynth_data()
-
-    def restore_data(self):
-        """Restore original data after resynthesis."""
-        self.fitness.restore_data()
+    def __repr__(self):
+        return "FitProblem(name=%s)" % self.name
 
     def valid(self, pvec):
         """Return true if the point is in the feasible region"""
-        return all(v in p.bounds for p, v in zip(self._parameters, pvec))
+        return all(v in p.prior for p, v in zip(self._parameters, pvec))
 
     def setp(self, pvec):
         """
@@ -324,7 +335,7 @@ def setp(self, pvec):
         # WARNING: don't try to conditionally update the model
         # depending on whether any model parameters have changed.
         # For one thing, the model_update below probably calls
-        # the subclass MultiFitProblem.model_update, which signals
+        # the subclass FitProblem.model_update, which signals
         # the individual models.  Furthermore, some parameters may
         # related to others via expressions, and so a dependency
         # tree needs to be generated.  Whether this is better than
@@ -344,8 +355,8 @@ def getp(self):
         return np.array([p.value for p in self._parameters], 'd')
 
     def bounds(self):
-        """Return the bounds fore each parameter a 2 x N array"""
-        limits = [p.bounds.limits for p in self._parameters]
+        """Return the bounds for each parameter as a 2 x N array"""
+        limits = [p.prior.limits for p in self._parameters]
         return np.array(limits, 'd').T if limits else np.empty((2, 0))
 
     def randomize(self, n=None):
@@ -365,73 +376,68 @@ def randomize(self, n=None):
             self.setp(self.randomize(n=1)[0])
             return   # Not returning anything since no n is requested
 
+        # TODO: apply hard limits on parameters
         target = self.getp()
         target[~np.isfinite(target)] = 1.
-        pop = [p.bounds.random(n, target=v)
+        pop = [p.prior.random(n, target=v)
                for p, v in zip(self._parameters, target)]
         return np.array(pop).T
 
-    def parameter_nllf(self):
-        """
-        Returns negative log likelihood of seeing parameters p.
-        """
-        s = sum(p.nllf() for p in self.bounded)
-        # print "; ".join("%s %g %g"%(p,p.value,p.nllf()) for p in
-        # self.bounded)
-        return s
-
-    def constraints_nllf(self):
-        """
-        Returns the cost of all constraints.
+    def nllf(self, pvec=None) -> float:
         """
-        return self.constraints() if self.constraints else 0.
+        compute the cost function for a new parameter set p.
 
-    def parameter_residuals(self):
-        """
-        Returns negative log likelihood of seeing parameters p.
-        """
-        return [p.residual() for p in self.bounded]
+        this is not simply the sum-squared residuals, but instead is the
+        negative log likelihood of seeing the data given the model parameters
+        plus the negative log likelihood of seeing the model parameters.  the
+        value is used for a likelihood ratio test so normalization constants
+        can be ignored.  there is an additional penalty value provided by
+        the model which can be used to implement inequality constraints.  any
+        penalty should be large enough that it is effectively excluded from
+        the parameter space returned from uncertainty analysis.
 
-    @property
-    def has_residuals(self):
+        the model is not actually calculated if any of the parameters are 
+        out of bounds, or any of the constraints are not satisfied, but
+        instead are assigned a value of *penalty_nllf*.  
+        this will prevent expensive models from spending time computing 
+        values in the unfeasible region.
         """
-        True if the underlying fitness function defines residuals.
-        """
-        return hasattr(self.fitness, 'residuals')
+        if pvec is not None:
+            if self.valid(pvec):
+                self.setp(pvec)
+            else:
+                return inf
 
-    def residuals(self):
-        r"""
-        Return the model residuals.
+        pparameter, pconstraints, pmodel, failing_constraints = self._nllf_components()
+        cost = pparameter + pconstraints + pmodel
+        # print(pvec, "cost=",pparameter,"+",pconstraints,"+",pmodel,"=",cost)
+        if isnan(cost):
+            # TODO: make sure errors get back to the user
+            # print "point evaluates to nan"
+            # print parameter.summarize(self._parameters)
+            return inf
+        return cost
 
-        If the model is defined by $y = f(x) + \epsilon$ for normally
-        distributed error in the measurement $y$ equal to
-        $\epsilon \sim N(0, \sigma^2)$, then residuals will be defined by
-        $R = (y - f(x))/\sigma$.  If the measurement uncertainty is not normal,
-        then the normal equivalent residuals should be defined so that the
-        Levenberg-Marquardt fit behaves reasonably, and the plot of
-        residuals gives an indication of which points are driving the fit.
+    def parameter_nllf(self) -> util.Tuple[float, util.List[str]]:
         """
-        if not hasattr(self.fitness, 'residuals'):
-            ## Fake residuals by using single nllf value as residual
-            #return np.asarray([np.sqrt(self.nllf())])
-            raise NotImplementedError("model does not define residuals")
-        return self.fitness.residuals()
-
-    def chisq(self):
+        Returns negative log likelihood of seeing parameters p.
         """
-        Return sum squared residuals normalized by the degrees of freedom.
-
-        In the context of a composite fit, the reduced chisq on the individual
-        models only considers the points and the fitted parameters within
-        the individual model.
+        failing = []
+        nllf = 0.0
+        for p in self._bounded:
+            p_nllf = p.nllf()
+            nllf += p_nllf
+            if p_nllf == np.inf:
+                failing.append(str(p))
+        # print "; ".join("%s %g %g"%(p,p.value,p.nllf()) for p in
+        # self._bounded)
+        return nllf, failing
 
-        Note that this does not include cost factors due to constraints on
-        the parameters, such as sample_offset ~ N(0,0.01).
+    def parameter_residuals(self):
+        """
+        Returns negative log likelihood of seeing parameters p.
         """
-        if hasattr(self.fitness, 'chisq'):
-            return self.fitness.chisq()
-        return np.sum(self.residuals() ** 2) / self.dof
-        # return 2*self.nllf()/self.dof
+        return [p.residual() for p in self._bounded]
 
     def chisq_str(self):
         """
@@ -445,72 +451,41 @@ def chisq_str(self):
         cost of the prior parameters and the cost of the penalty constraints
         in the total nllf.  The constraint value is displayed separately.
         """
-        pparameter, pconstraints, pmodel = self._nllf_components()
+        pparameter, pconstraints, pmodel, failing_constraints = self._nllf_components()
         chisq_norm, chisq_err = nllf_scale(self)
         chisq = pmodel * chisq_norm
         text = format_uncertainty(chisq, chisq_err)
         constraints = pparameter + pconstraints
         if constraints > 0.:
             text += " constraints=%g" % constraints
+        if len(failing_constraints) > 0:
+            text += " failing_constraints=%s" % str(failing_constraints)
 
         return text
 
-    def nllf(self, pvec=None):
-        """
-        compute the cost function for a new parameter set p.
-
-        this is not simply the sum-squared residuals, but instead is the
-        negative log likelihood of seeing the data given the model parameters
-        plus the negative log likelihood of seeing the model parameters.  the
-        value is used for a likelihood ratio test so normalization constants
-        can be ignored.  there is an additional penalty value provided by
-        the model which can be used to implement inequality constraints.  any
-        penalty should be large enough that it is effectively excluded from
-        the parameter space returned from uncertainty analysis.
-
-        the model is not actually calculated if the parameter nllf plus the
-        constraint nllf are bigger than *soft_limit*, but instead it is
-        assigned a value of *penalty_nllf*.  this will prevent expensive
-        models from spending time computing values in the unfeasible region.
-        """
-        if pvec is not None:
-            if self.valid(pvec):
-                self.setp(pvec)
-            else:
-                return inf
-
-        pparameter, pconstraints, pmodel = self._nllf_components()
-        cost = pparameter + pconstraints + pmodel
-        # print(pvec, "cost=",pparameter,"+",pconstraints,"+",pmodel,"=",cost)
-        if isnan(cost):
-            # TODO: make sure errors get back to the user
-            # print "point evaluates to nan"
-            # print parameter.summarize(self._parameters)
-            return inf
-        return cost
-
-    def _nllf_components(self):
+    def _nllf_components(self) -> util.Tuple[float, float, float, util.List[str]]:
         try:
-            pparameter = self.parameter_nllf()
+            pparameter, failing_parameter_constraints = self.parameter_nllf()
             if isnan(pparameter):
                 # TODO: make sure errors get back to the user
                 info = ["Parameter nllf is wrong"]
-                info += ["%s %g"%(p, p.nllf()) for p in self.bounded]
+                info += ["%s %g"%(p, p.nllf()) for p in self._bounded]
                 logging.error("\n  ".join(info))
-            pconstraints = self.constraints_nllf()
+            pconstraints, failing_constraints = self.constraints_nllf()
+            failing_constraints += failing_parameter_constraints
             # Note: for hard constraints (which return inf) avoid computing
             # model even if soft_limit is inf by using strict comparison
             # since inf <= inf is True but inf < inf is False.
-            pmodel = (self.model_nllf()
-                      if pparameter + pconstraints < self.soft_limit
-                      else self.penalty_nllf)
-            return pparameter, pconstraints, pmodel
+            penalty_nllf = self.penalty_nllf if self.penalty_nllf is not None else np.inf
+            pmodel = self.model_nllf() if len(failing_constraints) == 0 else penalty_nllf
+                     
+            return pparameter, pconstraints, pmodel, failing_constraints
         except Exception:
             # TODO: make sure errors get back to the user
             info = (traceback.format_exc(),
                     parameter.summarize(self._parameters))
             logging.error("\n".join(info))
-            return NaN, NaN, NaN
+            return NaN, NaN, NaN, []
 
     def __call__(self, pvec=None):
         """
@@ -523,157 +498,92 @@ def __call__(self, pvec=None):
         """
         return 2 * self.nllf(pvec) / self.dof
 
-    def show(self, _subs={}):
-        """Print the available parameters to the console as a tree."""
-        print(parameter.format(self.model_parameters(), freevars=_subs))
-        print("[chisq=%s, nllf=%g]" % (self.chisq_str(), self.nllf()))
-        #print(self.summarize())
-
     def summarize(self):
         """Return a table of current parameter values with range bars."""
         return parameter.summarize(self._parameters)
 
-    def labels(self):
+    def labels(self) -> util.List[str]:
         """Return the list of labels, one per fitted parameter."""
         return [p.name for p in self._parameters]
 
-    def save(self, basename):
-        """
-        Save the problem state for the current parameter set.
-
-        The underlying Fitness object *save* method is called, if it exists,
-        so that theory values can be saved in a format suitable to the problem.
-
-        Uses *basename* as the base of any files that are created.
-        """
-        if hasattr(self.fitness, 'save'):
-            self.fitness.save(basename)
-
-    def plot(self, p=None, fignum=None, figfile=None, view=None):
-        """
-        Plot the problem state for the current parameter set.
-
-        The underlying Fitness object *plot* method is called with *view*.
-        It should produce its plot on the current matplotlib figure.  This
-        method will add chisq to the plot and save it to a file.
-        """
-        if not hasattr(self.fitness, 'plot'):
-            return
-
-        import pylab
-        if fignum is not None:
-            pylab.figure(fignum)
-        if p is not None:
-            self.setp(p)
-        self.fitness.plot(view=view)
-        pylab.text(0.01, 0.01, 'chisq=%s' % self.chisq_str(),
-                   transform=pylab.gca().transAxes)
-        if figfile is not None:
-            pylab.savefig(figfile + "-model.png", format='png')
-
-    def cov(self):
-        """
-        Return the covariance matrix as computed by numdifftools from the
-        Hessian matrix for the problem at the current parameter values.
+    # noinspection PyAttributeOutsideInit
+    def model_reset(self):
         """
-        # TODO: remove from model
-        warnings.warn("use cov and stderr from FitDriver, not problem.",
-                      DeprecationWarning)
+        Prepare for the fit.
 
-        from . import lsqerror
-        H = lsqerror.hessian(self)
-        H, L = lsqerror.perturbed_hessian(H)
-        return lsqerror.chol_cov(L)
+        This sets the parameters and the bounds properties that the
+        solver is expecting from the fittable object.  We also compute
+        the degrees of freedom so that we can return a normalized fit
+        likelihood.
 
-    def stderr(self):
-        """
-        Return the 1-sigma uncertainty estimate for each parameter and the
-        correlation matrix *R* as computed from the covariance returned by
-        *cov*.
+        If the set of fit parameters changes, then model_reset must
+        be called.
         """
-        # TODO: remove from model
-        warnings.warn("use cov and stderr from FitDriver, not problem.",
-                      DeprecationWarning)
-
-        from . import lsqerror
-        c = self.cov()
-        return lsqerror.stderr(c), lsqerror.corr(c)
-
-    def __getstate__(self):
-        return (self.fitness, self.partial, self.name, self.penalty_nllf,
-                self.soft_limit, self.constraints)
-
-    def __setstate__(self, state):
-        self.fitness, self.partial, self.name, self.penalty_nllf, \
-            self.soft_limit, self.constraints = state
-        self.model_reset()
-
-class MultiFitProblem(BaseFitProblem):
-    """
-    Weighted fits for multiple models. See :func:`FitProblem` for an
-    explanation of weights.
-    """
-    def __init__(self, models, weights=None, name=None,
-                 constraints=None,
-                 soft_limit=np.inf, penalty_nllf=1e6,
-                 freevars=None):
-        self.partial = False
-        self.constraints = constraints
-        if freevars is None:
-            names = ["M%d" % i for i, _ in enumerate(models)]
-            freevars = parameter.FreeVariables(names=names)
-        self.freevars = freevars
-        self._models = [BaseFitProblem(m, partial=True) for m in models]
-        if weights is None:
-            weights = [1 for _ in models]
-        self.weights = weights
-        self.penalty_nllf = penalty_nllf
-        self.soft_limit = soft_limit
-        self.set_active_model(0)  # Set the active model to model 0
-        self.model_reset()
-        self.name = name
-
-    @property
-    def models(self):
-        """Iterate over models, with free parameters set from model values"""
-        for i, f in enumerate(self._models):
-            self.freevars.set_model(i)
-            yield f
-        # Restore the active model after cycling
-        self.freevars.set_model(self._active_model_index)
-
-    # noinspection PyAttributeOutsideInit
-    def set_active_model(self, i):
-        """Use free parameters from model *i*"""
-        self._active_model_index = i
-        self.active_model = self._models[i]
-        self.freevars.set_model(i)
-
-    def model_parameters(self):
-        """Return parameters from all models"""
-        pars = {'models': [f.model_parameters() for f in self.models]}
-        free = self.freevars.parameters()
-        if free:
-            pars['freevars'] = free
-        return pars
-
-    def to_dict(self):
-        return {
-            'type': type(self).__name__,
-            'name': self.name,
-            'models': to_dict(self._models),
-            'weights': self.weights,
-            'partial': self.partial,
-            'soft_limit': self.soft_limit,
-            'penalty_nllf': self.penalty_nllf,
-            # TODO: constraints may be a function.
-            'constraints': to_dict(self.constraints),
-            'freevars': to_dict(self.freevars),
-        }
+        # print self.model_parameters()
+        all_parameters = parameter.unique(self.model_parameters())
+        # print "all_parameters",all_parameters
+        # for p in all_parameters:
+        #     if hasattr(p, 'reset_prior'):
+        #         p.reset_prior()  # no constraints
+        #     else:
+        #         raise ValueError(f"{p} does not have prior")
+        broken = []
+        for p in all_parameters:
+            # slot = p.slot
+            # value = p.value
+            
+            # TODO: this is a shim to accomodate Expression, Calculation etc. being
+            # put into attributes that have type "Parameter" (in user scripts)
+            # Do we cause those scripts to break instead?  
+            # Or do we autoconvert to .equals()?
+            if hasattr(p, 'add_prior'):
+                p.add_prior(p.distribution, bounds=p.bounds, limits=p.limits)
+
+            # TODO: currently this logic for showing breaking constraints to the user 
+            # is in the chisq_str execution path instead of on model_reset - should it be here instead?
+
+            # While we are walking all parameters check which constraints aren't satisfied
+            # Build up a list of strings to help the user initialize the model correctly
+            if hasattr(p, 'limits'):
+                value = p.value
+                if (p.limits[0] > value) or (value > p.limits[1]):
+                    broken.append(f"{p}={value} is outside {p.limits}")
+                elif not isfinite(p.prior.nllf(value)):
+                    broken.append(f"{p}={value} is outside {p.prior}")
+
+        broken.extend([f"constraint {c} is unsatisfied" for c in self.constraints if float(c) == inf])
+        if self._constraints_function() == inf:
+             broken.append("user constraint function is unsatisfied")
+        if len(broken) > 0:
+            warnings.warn("Unsatisfied constraints: [%s]" % (",\n".join(broken)))
+        self.broken_constraints = broken
+
+        # TODO: shimmed to allow non-Parameter in Parameter attribute spots.
+        pars = []
+        pars_by_id = {}
+        for p in all_parameters:
+            if hasattr(p, 'id'):
+                pars_by_id[p.id] = p
+            slot = getattr(p, 'slot', None)
+            if isinstance(slot, Variable) and not p.fixed:
+                pars.append(p)
+        self._parameters = pars
+        self._parameters_by_id = pars_by_id
+        # TODO: shimmed to allow non-Parameter in Parameter attribute spots.
+        self._bounded = [p for p in all_parameters if hasattr(p, 'has_prior') and p.has_prior()]
+        self.dof = self.model_points()
+        self.dof -= len(self._parameters)
+        if self.dof <= 0:
+            raise ValueError("Need more data points than fitting parameters")
+        #self.constraints = pars.constraints()
+        # Find the constraints on variables and expressions that we need to compute
+        # parameter_constraints = [p.slot for p in all_parameters if isinstance(p.slot, Variable) and p.has_prior()]
+        # expression_constraints = [p.slot for p in all_parameters if isinstance(p.slot, Expression) and p.has_prior()]
+        # self._all_constraints = parameter_constraints + expression_constraints
 
     def model_points(self):
         """Return number of points in all models"""
-        return sum(f.model_points() for f in self.models)
+        return sum(f.numpoints() for f in self.models)
 
     def model_update(self):
         """Let all models know they need to be recalculated"""
@@ -685,16 +595,28 @@ def model_update(self):
         # computation is cheap before jumping to a more expensive
         # subspace.  SrFit does this.
         for f in self.models:
-            f.model_update()
+            if hasattr(f, 'update'):
+                f.update()
 
     def model_nllf(self):
         """Return cost function for all data sets"""
-        return sum(w**2*f.model_nllf() for w, f in zip(self.weights, self.models))
+        return sum(w**2*f.nllf() for w, f in zip(self.weights, self.models))
 
-    def constraints_nllf(self):
+    def constraints_nllf(self) -> util.Tuple[float, util.List[str]]:
         """Return the cost function for all constraints"""
-        return (sum(f.constraints_nllf() for f in self.models)
-                + BaseFitProblem.constraints_nllf(self))
+        failing = []
+        nllf = 0.0
+        nllf = self._constraints_function()
+        if nllf == np.inf:
+            failing.append("user constraints function")
+        for c in self.constraints:
+            # TODO: convert to list of residuals for Levenberg-Marquardt
+            c_nllf = float(c)**2
+            nllf += c_nllf
+            if c_nllf > 0:
+                failing.append(str(c))
+
+        return nllf, failing
 
     def simulate_data(self, noise=None):
         """Simulate data with added noise"""
@@ -716,7 +638,7 @@ def has_residuals(self):
         """
         True if all underlying fitness functions define residuals.
         """
-        return all(f.has_residuals for f in self.models)
+        return all(hasattr(f, 'residuals') for f in self.models)
 
     def residuals(self):
         resid = np.hstack([w * f.residuals()
@@ -725,13 +647,14 @@ def residuals(self):
 
     def save(self, basename):
         for i, f in enumerate(self.models):
-            f.save(basename + "-%d" % (i + 1))
+            if hasattr(f, 'save'):
+                f.save(basename + "-%d" % (i + 1))
 
     def show(self):
         for i, f in enumerate(self.models):
-            print("-- Model %d %s" % (i, f.name))
+            print("-- Model %d %s" % (i, getattr(f, 'name', '')))
             subs = self.freevars.get_model(i) if self.freevars else {}
-            f.show(_subs=subs)
+            show_parameters(f, subs=subs)
         print("[overall chisq=%s, nllf=%g]" % (self.chisq_str(), self.nllf()))
 
     def plot(self, p=None, fignum=1, figfile=None, view=None):
@@ -739,8 +662,14 @@ def plot(self, p=None, fignum=1, figfile=None, view=None):
         if p is not None:
             self.setp(p)
         for i, f in enumerate(self.models):
-            f.plot(fignum=i + fignum, view=view)
+            if not hasattr(f, 'plot'):
+                continue
+            f.plot(view=view)
+            pylab.figure(i + fignum)
+            f.plot(view=view)
             pylab.suptitle('Model %d - %s' % (i, f.name))
+            pylab.text(0.01, 0.01, 'chisq=%s' % chisq_str(f),
+                   transform=pylab.gca().transAxes)
             if figfile is not None:
                 pylab.savefig(figfile + "-model%d.png" % i, format='png')
 
@@ -755,7 +684,7 @@ def __setstate__(self, state):
 
 # TODO: consider adding nllf_scale to FitProblem.
 ONE_SIGMA = 0.68268949213708585
-def nllf_scale(problem):
+def nllf_scale(problem: FitProblem):
     r"""
     Return the scale factor for reporting the problem nllf as an approximate
     normalized chisq, along with an associated "uncertainty".  The uncertainty
@@ -780,7 +709,7 @@ def nllf_scale(problem):
         npars = max(len(problem.getp()), 1)
         return 2./dof, chi2.ppf(ONE_SIGMA, npars)/dof
 
-def load_problem(filename, options=None):
+def load_problem(filename, options=None) -> FitProblem:
     """
     Load a problem definition from a python script file.
 
diff --git a/bumps/fitters.py b/bumps/fitters.py
index 01a67a7e..12f14761 100644
--- a/bumps/fitters.py
+++ b/bumps/fitters.py
@@ -608,7 +608,7 @@ def _residuals(self, p, fjac=None):
         residuals = np.hstack(
             (self.problem.residuals().flat, self.problem.parameter_residuals()))
         # Tally costs for broken constraints
-        extra_cost = self.problem.constraints_nllf()
+        extra_cost, failing_constraints = self.problem.constraints_nllf()
         # Spread the cost over the residuals.  Since we are smoothly increasing
         # residuals as we leave the boundary, this should push us back into the
         # boundary (within tolerance) during the lm fit.
@@ -681,7 +681,8 @@ def _bounded_residuals(self, p):
         residuals = np.hstack(
             (self.problem.residuals().flat, self.problem.parameter_residuals()))
         # Tally costs for straying outside the boundaries plus other costs
-        extra_cost = stray_cost + self.problem.constraints_nllf()
+        constraints_cost, failing_constraints = self.problem.constraints_nllf()
+        extra_cost = stray_cost + constraints_cost
         # Spread the cost over the residuals.  Since we are smoothly increasing
         # residuals as we leave the boundary, this should push us back into the
         # boundary (within tolerance) during the lm fit.
diff --git a/bumps/gui/app_panel.py b/bumps/gui/app_panel.py
index 7f408a00..3455a4b7 100644
--- a/bumps/gui/app_panel.py
+++ b/bumps/gui/app_panel.py
@@ -624,7 +624,9 @@ def save_model(self, path):
         except ImportError:
             from pickle import dump
         try:
-            if hasattr(self.model, 'save_json'):
+            if hasattr(plugin, 'save_json'):
+                plugin.save_json(self.model, path)
+            elif hasattr(self.model, 'save_json'):
                 self.model.save_json(path)
             with open(path,'wb') as fid:
                 dump(self.model, fid)
@@ -711,6 +713,8 @@ def set_model(self, model):
         # Enable appropriate toolbar items.
         self.tb.EnableTool(self.tb_start.GetId(), True)
         self.tb.EnableTool(self.tb_stop.GetId(), True)
+        if hasattr(model, 'broken_constraints'):
+            signal.log_message(message="Unsatisfied constraints: \n\t%s" % (",\n\t".join(model.broken_constraints)))
         if hasattr(model, 'path'):
             signal.log_message(message="loaded "+model.path)
             self.GetTopLevelParent().SetTitle("Bumps: %s"%model.name)
diff --git a/bumps/gui/data_view.py b/bumps/gui/data_view.py
index 04a11a78..8e90d141 100755
--- a/bumps/gui/data_view.py
+++ b/bumps/gui/data_view.py
@@ -12,7 +12,7 @@
 # The Figure object is used to create backend-independent plot representations.
 from matplotlib.figure import Figure
 
-from ..fitproblem import MultiFitProblem
+from ..fitproblem import FitProblem
 
 from .util import EmbeddedPylab
 
@@ -157,7 +157,7 @@ def redraw(self, reset=False):
                     break
 
             # Preform the calculation
-            if isinstance(self.problem,MultiFitProblem):
+            if isinstance(self.problem,FitProblem):
                 #print "n=",len(self.problem.models)
                 for p in self.problem.models:
                     self._precalc(p)
@@ -179,10 +179,10 @@ def redraw(self, reset=False):
                 pylab.clf() # clear the canvas
                 #shift=20 if self.view == 'log' else 0
                 shift=0
-                if isinstance(self.problem, MultiFitProblem):
+                if isinstance(self.problem, FitProblem):
                     for i,p in enumerate(self.problem.models):
                         #if hasattr(p.fitness,'plot'):
-                        p.fitness.plot(view=self.view)
+                        p.plot(view=self.view)
                         if self._cancel_calculate: break
                     pylab.text(0.01, 0.01,
                                'chisq=%s' % self.problem.chisq_str(),
diff --git a/bumps/gui/parameter_view.py b/bumps/gui/parameter_view.py
index 80ccac4c..88b0a13d 100755
--- a/bumps/gui/parameter_view.py
+++ b/bumps/gui/parameter_view.py
@@ -31,7 +31,7 @@
 from uuid import uuid4 as uuid_generate
 import sys
 
-from ..parameter import BaseParameter
+from ..parameter import Parameter
 from .util import nice
 from . import signal
 
@@ -115,7 +115,7 @@ def IsContainer(self, item):
 
         # Return False if it is a leaf
         node = self.ItemToObject(item)
-        return not isinstance(node["value"], BaseParameter)
+        return not isinstance(node["value"], Parameter)
 
     def GetParent(self, item):
         # Return the item which is this item's parent.
@@ -161,14 +161,14 @@ def GetValue(self, item, col):
             return mapper[col]
 
 
-        elif isinstance(par, BaseParameter):
+        elif isinstance(par, Parameter):
             if par.fittable:
                 if par.fixed:
                     fitted = False
                     low, high = '', ''
                 else:
                     fitted = True
-                    low, high = (str(v) for v in par.bounds.limits)
+                    low, high = (str(v) for v in par.prior.limits)
             else:
                 fitted = False
                 low, high = '', ''
@@ -194,7 +194,7 @@ def SetValue(self, value, item, col):
 
         node = self.ItemToObject(item)
         par = node["value"]
-        if isinstance(par, BaseParameter):
+        if isinstance(par, Parameter):
             if col == 0:
                 if par.fittable:
                     par.fixed = not value
@@ -203,14 +203,14 @@ def SetValue(self, value, item, col):
             elif col == 3:
                 if value == '': return
                 low = float(value)
-                high = par.bounds.limits[1]
-                if low != par.bounds.limits[0]:
+                high = par.prior.limits[1]
+                if low != par.prior.limits[0]:
                     par.range(low, high)
             elif col == 4:
                 if value == '': return
                 high = float(value)
-                low = par.bounds.limits[0]
-                if high != par.bounds.limits[1]:
+                low = par.prior.limits[0]
+                if high != par.prior.limits[1]:
                     par.range(low, high)
 
         if col == 0:
@@ -314,14 +314,14 @@ def params_to_dict(params):
             ref[k] = params_to_dict(params[k])
     elif isinstance(params, tuple) or isinstance(params, list):
         ref = [params_to_dict(v) for v in params]
-    elif isinstance(params, BaseParameter):
+    elif isinstance(params, Parameter):
         if params.fittable:
             if params.fixed:
                 fitted = 'No'
                 low, high = '', ''
             else:
                 fitted = 'Yes'
-                low, high = (str(v) for v in params.bounds.limits)
+                low, high = (str(v) for v in params.prior.limits)
         else:
             fitted = ''
             low, high = '', ''
@@ -346,7 +346,7 @@ def params_to_list(params, parent_uuid=None, output=None, path='M', links=None):
             #output.append(new_item)
             new_id = None
             params_to_list(v, parent_uuid=new_id, output=output, path=path+"[%d]" % (i,), links=links)
-    elif isinstance(params, BaseParameter):
+    elif isinstance(params, Parameter):
         link_path = links.get(id(params), "")
         if link_path == "":
             links[id(params)] = path
diff --git a/bumps/gui/summary_view.py b/bumps/gui/summary_view.py
index 0d9e23c3..d31f1272 100755
--- a/bumps/gui/summary_view.py
+++ b/bumps/gui/summary_view.py
@@ -173,7 +173,7 @@ def __init__(self, parent, parameter, model):
         self.parameter = parameter
         self.model = model
 
-        self.low, self.high = (v for v in self.parameter.bounds.limits)
+        self.low, self.high = (v for v in self.parameter.prior.limits)
 
         text_hbox = wx.BoxSizer(wx.HORIZONTAL)
 
@@ -204,14 +204,14 @@ def __init__(self, parent, parameter, model):
         self.update_slider()
 
     def update_slider(self):
-        slider_pos = int(self.parameter.bounds.get01(self.parameter.value)*NUMTICKS)
+        slider_pos = int(self.parameter.prior.get01(self.parameter.value)*NUMTICKS)
         # Add line below if get01 doesn't protect against values out of range.
         #slider_pos = min(max(slider_pos,0),100)
         self.slider.SetValue(slider_pos)
         self.value.SetLabel(VALUE_FORMAT.format(nice(self.parameter.value, digits=VALUE_PRECISION)))
 
         # Update new min and max range of values if changed.
-        newlow, newhigh = (v for v in self.parameter.bounds.limits)
+        newlow, newhigh = (v for v in self.parameter.prior.limits)
         if newlow != self.low:
             self.min_range.SetLabel(VALUE_FORMAT.format(newlow))
 
@@ -220,7 +220,7 @@ def update_slider(self):
 
     def OnScroll(self, event):
         value = self.slider.GetValue()
-        new_value  = self.parameter.bounds.put01(value/NUMTICKS)
+        new_value  = self.parameter.prior.put01(value/NUMTICKS)
         nice_new_value = nice(new_value, digits=VALUE_PRECISION)
         self.parameter.value = nice_new_value
         self.value.SetLabel(VALUE_FORMAT.format(nice_new_value))
diff --git a/bumps/lsqerror.py b/bumps/lsqerror.py
index 9d37d8b6..f34f66cd 100644
--- a/bumps/lsqerror.py
+++ b/bumps/lsqerror.py
@@ -68,6 +68,12 @@ def jacobian(problem, p=None, step=None):
     point value p_j if the range on parameter j is infinite.
 
     The current point is preserved.
+
+    Note that the problem.residuals() method should not reuse memory for the
+    returned value otherwise the derivative calculation (f(x+dx) - f(x))/dx
+    will always be zero. The returned value need not be 1D, but it should be
+    contiguous so that it can be reshaped to 1D without an extra copy. This
+    will only be an issue for very large datasets.
     """
     p_init = problem.getp()
     if p is None:
@@ -76,7 +82,10 @@ def jacobian(problem, p=None, step=None):
     bounds = getattr(problem, 'bounds', lambda: None)()
     def f(p):
         problem.setp(p)
-        return problem.residuals()
+        # Return residuals as a vector even if f(x) returns a matrix otherwise
+        # we cannot build a stacked Jacobian. We use reshape() rather than
+        # flatten since this will avoid an unnecessary copy.
+        return np.reshape(problem.residuals(), -1)
     J = _jacobian_forward(f, p, bounds, eps=step)
     #J = nd.Jacobian(problem.residuals)(p)
     problem.setp(p_init)
@@ -86,7 +95,6 @@ def _jacobian_forward(f, p, bounds, eps=None):
     n = len(p)
     # TODO: default to double precision epsilon
     step = 1e-4 if eps is None else np.sqrt(eps)
-    fx = f(p)
 
     #print("p",p,"step",step)
     h = abs(p)*step
@@ -95,9 +103,11 @@ def _jacobian_forward(f, p, bounds, eps=None):
         h[h+p > bounds[1]] *= -1.0  # step backward if forward step is out of bounds
     ee = np.diag(h)
 
+    fx = f(p) # Maybe fx.copy() to protect against reuse
     J = []
     for i in range(n):
-        J.append((f(p + ee[i, :]) - fx)/h[i])
+        fx_plus = f(p + ee[i, :])
+        J.append((fx_plus - fx)/h[i])
     return np.vstack(J).T
 
 def _jacobian_central(f, p, bounds, eps=None):
@@ -114,7 +124,9 @@ def _jacobian_central(f, p, bounds, eps=None):
 
     J = []
     for i in range(n):
-        J.append((f(p + ee[i, :]) - f(p - ee[i, :])) / (2.0*h[i]))
+        fx_minus = f(p - ee[i, :]) # Maybe fx.copy() to protect against reuse
+        fx_plus = f(p + ee[i, :])
+        J.append((fx_plus - fx_minus) / (2.0*h[i]))
     return np.vstack(J).T
 
 
diff --git a/bumps/mpfit.py b/bumps/mpfit.py
index 736547cc..0b346f8a 100644
--- a/bumps/mpfit.py
+++ b/bumps/mpfit.py
@@ -141,7 +141,7 @@ def myfunct(p, fjac=None, x=None, y=None, err=None)
     # stop the calculation.
     status = 0
     if (dojac):
-       pderiv = numpy.zeros([len(x), len(p)], numpy.float)
+       pderiv = numpy.zeros([len(x), len(p)], numpy.float64)
        for j in range(len(p)):
          pderiv[:,j] = FGRAD(x, p, j)
     else:
@@ -279,7 +279,7 @@ def myfunct(p, fjac=None, x=None, y=None, err=None)
 
    import mpfit
    import numpy
-   x = numpy.arange(100, numpy.float)
+   x = numpy.arange(100, numpy.float64)
    p0 = [5.7, 2.2, 500., 1.5, 2000.]
    y = ( p[0] + p[1]*[x] + p[2]*[x**2] + p[3]*numpy.sqrt(x) +
          p[4]*numpy.log(x))
@@ -897,8 +897,8 @@ def iterfunct(myfunct, p, iter, fnorm, functkw=None,
             self.errmsg = 'ERROR: either P or PARINFO(*)["value"] must be supplied.'
             return
 
-      ## Make sure parameters are numpy arrays of type numpy.float
-      xall = numpy.asarray(xall, numpy.float)
+      ## Make sure parameters are numpy arrays of type numpy.float64
+      xall = numpy.asarray(xall, numpy.float64)
 
       npar = len(xall)
       self.fnorm  = -1.
@@ -995,8 +995,8 @@ def iterfunct(myfunct, p, iter, fnorm, functkw=None,
          if (len(wh) > 0): return
          self.errmsg = ''
 
-      # Make sure x is a numpy array of type numpy.float
-      x = numpy.asarray(x, numpy.float)
+      # Make sure x is a numpy array of type numpy.float64
+      x = numpy.asarray(x, numpy.float64)
 
       [self.status, fvec] = self.call(fcn, self.params, functkw)
       if (self.status < 0):
@@ -1329,14 +1329,14 @@ def iterfunct(myfunct, p, iter, fnorm, functkw=None,
 
             ## Fill in actual covariance matrix, accounting for fixed
             ## parameters.
-            self.covar = numpy.zeros([nn, nn], numpy.float)
+            self.covar = numpy.zeros([nn, nn], numpy.float64)
             for i in range(n):
                indices = ifree+ifree[i]*n
                numpy.put(self.covar, indices, cv[:,i])
 
             ## Compute errors in parameters
             catch_msg = 'computing parameter errors'
-            self.perror = numpy.zeros(nn, numpy.float)
+            self.perror = numpy.zeros(nn, numpy.float64)
             d = numpy.diagonal(self.covar).copy()
             wh, = numpy.nonzero(d >= 0)
             if len(wh) > 0:
@@ -1408,9 +1408,9 @@ def parinfo(self, parinfo=None, key='a', default=None, n=0):
       test = default
       if isinstance(default, list): test=default[0]
       if isinstance(test, int):
-         values = numpy.asarray(values, numpy.int)
+         values = numpy.asarray(values, numpy.int32)
       elif isinstance(test, float):
-         values = numpy.asarray(values, numpy.float)
+         values = numpy.asarray(values, numpy.float64)
       return(values)
 
 
@@ -1482,7 +1482,7 @@ def fdjac2(self, fcn, x, fvec, step=None, ulimited=None, ulimit=None, dside=None
       ## Compute analytical derivative if requested
       if (autoderivative == 0):
          mperr = 0
-         fjac = numpy.zeros(nall, numpy.float)
+         fjac = numpy.zeros(nall, numpy.float64)
          numpy.put(fjac, ifree, 1.0)  ## Specify which parameters need derivatives
          [status, fp, pderiv] = self.call(fcn, xall, functkw, fjac=fjac)
 
@@ -1502,7 +1502,7 @@ def fdjac2(self, fcn, x, fvec, step=None, ulimited=None, ulimit=None, dside=None
             fjac.shape = [m, n]
             return(fjac)
 
-      fjac = numpy.zeros([m, n], numpy.float)
+      fjac = numpy.zeros([m, n], numpy.float64)
 
       h = eps * abs(x)
 
@@ -1687,7 +1687,7 @@ def qrfac(self, a, pivot=0):
       n = sz[1]
 
       ## Compute the initial column norms and initialize arrays
-      acnorm = numpy.zeros(n, numpy.float)
+      acnorm = numpy.zeros(n, numpy.float64)
       for j in range(n):
          acnorm[j] = self.enorm(a[:,j])
       rdiag = acnorm.copy()
diff --git a/bumps/names.py b/bumps/names.py
index c8692785..f29391c1 100644
--- a/bumps/names.py
+++ b/bumps/names.py
@@ -22,16 +22,12 @@
 - :class:`Curve <bumps.curve.Curve>` for defining models from functions
 - :class:`PoissonCurve <bumps.curve.PoissonCurve>` for modelling data with Poisson uncertainty
 - :class:`PDF <bumps.pdfwrapper.PDF>` for fitting a probability distribution directly
-- :func:`FitProblem <bumps.fitproblem.FitProblem>` for defining the fit (see
-    :class:`BaseFitProblem <bumps.fitproblem.BaseFitProblem>` or
-    :class:`MultiFitProblem <bumps.fitproblem.MultiFitProblem>` for details,
-    depending on whether you are fitting a single model or multiple models
-    simultaneously).
+- :class:`FitProblem <bumps.fitproblem.FitProblem>` for defining the fit.
 """
 
 #__all__ = [ 'sys', 'np', 'inf', 'pmath',
 #    'Parameter', 'FreeVariables', 'Distribution', 'PDF', 'Curve', 'PoissonCurve',
-#        'FitProblem', 'MultiFitProblem' ]
+#        'FitProblem' ]
 
 import sys
 import numpy as np
@@ -48,6 +44,6 @@
 from .bounds import Distribution
 from .pdfwrapper import PDF, VectorPDF, DirectProblem
 from .curve import Curve, PoissonCurve
-from .fitproblem import FitProblem, MultiFitProblem
+from .fitproblem import FitProblem, Fitness
 from .fitters import fit
 from .util import relative_import
diff --git a/bumps/parameter.py b/bumps/parameter.py
index 75bff124..c09e7c87 100644
--- a/bumps/parameter.py
+++ b/bumps/parameter.py
@@ -13,16 +13,29 @@
 #__all__ = [ 'Parameter']
 import operator
 import sys
-from six.moves import reduce
+import builtins
+from functools import reduce
 import warnings
 from copy import copy
-import math
+import uuid
 from functools import wraps
+from enum import Enum
+
+from typing import (
+    Type, TypeVar, Optional, Any, Union, Dict, Callable,
+    Tuple, List, Sequence)
+from .util import Literal
 
 import numpy as np
 from numpy import inf, isinf, isfinite
 
 from . import bounds as mbounds
+from . import pmath
+from .util import field, field_desc, schema, dataclass
+
+BoundsType = mbounds.BoundsType
+
+ValueType = Union['Expression', 'Parameter', 'Calculation', float]
 
 # TODO: avoid evaluation of subexpressions if parameters do not change.
 # This is especially important if the subexpression invokes an expensive
@@ -65,23 +78,252 @@ def to_dict(p):
         return str(p)
 
 
-class BaseParameter(object):
+@schema(init=True)
+class Uniform:
+    """ Uniform distribution with hard boundaries """
+
+
+@schema()
+class Normal:
+    """ Normal distribution (Gaussian) """
+    std: float = field_desc("standard deviation (1-sigma)")
+    mean: float = field_desc("center of the distribution")
+    def __init__(self, std: float, mean: float):
+        self.std = std
+        self.mean = mean
+
+# Leave out of schema for now.
+# TODO: determine if this is used by anyone
+#@schema(init=True)
+class UniformSoftBounded:
+    """ Uniform distribution with error-function PDF on boundaries """
+    std: float = field_desc("width of the edge distribution")
+
+
+DistributionType = Union[Uniform, Normal] #, UniformSoftBounded]
+
+class OperatorMixin:
     """
-    Root of the parameter class, defining arithmetic on parameters
+    The set of operations that can be performed on parameter-like objects
+    Parameter, Constant, Expression.
+
+    These include: +, -, *, /, //, **, abs, float, int
+
+    Also, numpy math functions: sin, cos, tan, ...
+
+    Much like abs(obj) => obj.__abs__(), np.sin(obj) => obj.sin()
     """
+    # float(value) is special: it returns the current value rather than
+    # becoming part of the parameter expression.
+    value: float
+    def __float__(self):
+        return float(self.value)
+    def __int__(self):
+        return int(self.value)
+    def __bool__(self):
+        # Note: __bool__ must return true or false, so we can't handle
+        # lazy constraint expressions like not a, a or b, a and b.
+        raise TypeError("use (p != 0) to test against zero")
+    ...  # operators and functions will be filled in later
 
-    # Parameters are fixed unless told otherwise
-    fixed = True
-    fittable = False
-    discrete = False
-    _bounds = mbounds.Unbounded()
-    name = None
-    value = None # value is an attribute of the derived class
+class ValueProtocol(OperatorMixin):
+    """
+    Values can be combined to form expressions
+    Provide a suite of operators for creating parameter expressions.
+    """
+    fittable: bool = False
+    fixed: bool = True
+    value: float
+    # TODO: Do values have names? Or do the names belong to the model parameter?
+    #name: str
+    # TODO: are priors on the parameter or on the value?
+    #bounds: Optional[BoundsType] = None
+    def parameters(self) -> List["Parameter"]:
+        # default implementation:
+        return []
+
+
+@schema()
+class Calculation(ValueProtocol): # the name Function is taken (though deprecated)
+    """
+    A Parameter with a model-specific, calculated value.
+    The function used to calculate this value should be well-documented in the 
+    description field, e.g.
+    Stack.thickness: description = "a sum of the thicknesses of all layers in the stack"
+
+    """
+    description: str
+    _function: Callable[[], float]  # added by the model; not serialized
+    def __init__(self, description: str=""):
+        self.description = description
+    @property
+    def value(self):
+        return self._function()
+    def __float__(self):
+        return self.value
+    def set_function(self, function):
+        self._function = function
+
+
+class SupportsPrior:
+    prior: Optional[BoundsType]
+    limits: Tuple[float, float]
+    distribution: DistributionType
+    bounds: BoundsType
+
+    def reset_prior(self):
+        self.prior = None
+
+    def has_prior(self):
+        return (
+            self.prior is not None
+            and not isinstance(self.prior, mbounds.Unbounded)
+            and self.prior.limits != (-np.inf, np.inf))
+    
+    def add_prior(self, distribution: Optional[DistributionType]=None, bounds: Optional[BoundsType]=None, limits: Optional[Tuple[float, float]]=None):
+        # use self values if they are found:
+        if distribution is None and self.distribution is not None:
+            distribution = self.distribution
+        if bounds is None and self.bounds is not None:
+            bounds = self.bounds
+        if limits is None:
+            if self.limits is not None:
+                limits = self.limits
+            else:
+                limits = (-inf, inf)
+
+        if bounds is not None:
+            # get the intersection of the limits here.
+            limits = (
+                np.clip(limits[0], *bounds),
+                np.clip(limits[1], *bounds)
+            )
+        if isinstance(distribution, Normal):
+            if limits == (-inf, inf):
+                prior = mbounds.BoundedNormal(mean=distribution.mean, std=distribution.std, limits=limits)
+            else:
+                prior = mbounds.Normal(mean=distribution.mean, std=distribution.std)
+        elif isinstance(distribution, UniformSoftBounded):
+            lo, hi = limits
+            prior = mbounds.SoftBounded(lo=lo, hi=hi, std=distribution.std)
+        elif isinstance(distribution, Uniform):
+            lo, hi = limits
+            if isinf(lo) and isinf(hi):
+                prior = mbounds.Unbounded()
+            elif isinf(lo):
+                prior = mbounds.BoundedAbove(hi)
+            elif isinf(hi):
+                prior = mbounds.BoundedBelow(lo)
+            else:
+                prior = mbounds.Bounded(lo, hi)
+        else:
+            raise ValueError("no distribution found matching %s" % (str(distribution)))
+        
+        self.prior = prior
+
+
+@schema(classname="Parameter", eq=False)
+class ParameterSchema:
+    """
+    Saved state for a slot in the parameter model.
+    """
+    id: str = field(metadata={"format": "uuid"})
+    name: Optional[str] = field(default=None, init=False)
+    fixed: bool = True
+    slot: Union['Variable', ValueType]
+    limits: Tuple[Union[float, Literal["-inf"]], Union[float, Literal["inf"]]] = (-inf, inf)
+    bounds: Optional[Tuple[Union[float, Literal["-inf"]], Union[float, Literal["inf"]]]] = None
+    distribution: DistributionType = field(default_factory=Uniform)
+    discrete: bool = False
+
+class Parameter(ValueProtocol, ParameterSchema, SupportsPrior):
+    """
+    A parameter is a container for a symbolic value.
+
+    Parameters have a prior probability, as set by a bounds constraint:
+
+        import numpy as np
+        from scipy.stats.distributions import lognorm
+        from bumps.parameter import Parameter
+
+        p = Parameter(3)
+        p.pmp(10)               # 3 +/- 10% uniform
+        p.pmp(-5,10)            # 3 in [2.85, 3.30] uniform
+        p.pm(2)                 # 3 +/- 2 uniform
+        p.pm(-1,2)              # 3 in [2,5] uniform
+        p.range(0,5)            # 3 in [0,5] uniform
+        p.dev(2)                # 3 +/- 2 gaussian
+        p.soft_range(2,5,2)     # 3 in [2,5] uniform with gauss wings
+        p.dev(2, limits=(0,6))  # 3 +/- 2 truncated gaussian
+        p.pdf(lognorm(3, 1))    # lognormal centered on 3, width 1.
+
+    Parameters have hard limits on the possible values, dictated by the model.
+    These bounds apply in addition to any other bounds.
+
+    Parameters can be constrained to be equal to another parameter or
+    parameter expression:
+
+        a, b = Parameter(3), Parameter(4)
+        p = Parameter(limits=(6, 10))
+        p.equals(a+b)
+        assert p.nllf() == 0.   # within the bounds
+        a.value = 20
+        assert np.isinf(p.nllf()) # out of bounds
+
+    Constraints on the computed value follow from the constraints on the
+    underlying parameters in addition to any hard limits on the parameter
+    value given by the model.
+
+    **Inputs**
+
+    *value* can be a constant, a variable, an expression or a link to
+    another parameter.
+
+    *bounds* are user-supplied limits on the parameter value within the model.
+    If bounds are supplied then the parameter defaults to fittable.
+
+    *distribution* is one of Uniform, Normal or UniformSoftBounded classes
+
+    *fixed* is True if the parameter is fixed, even if bounds are supplied.
+
+    *name* is the label associated with the parameter in plots. The names
+    need not be unique, but it will be confusing if there are duplicates.
+    The name will usually correspond to the role of the parameter in the
+    model. For models with sequences (e.g., layer numbers), try using a
+    layer name (e.g., based on the material in the layer) rather than a layer
+    number for parameters in that layer. This will make it easier for the
+    user to associate the parameters displayed at the end of the the fit
+    with the layer in the model. Also, when exploring the space of models,
+    the parameter names will be preserved even if a new layer is introduced
+    before the existing layers. That will allow the parameters from the
+    previous fit to be easily used as a seed for the fit to the new model.
+
+    *id* must be a unique identifier associated with the parameter. This
+    is used to link parameters on save and reload.
+
+    *limits* are hard limits on the parameter value within the model. Separate
+    from the prior distribution on a random variable provided by the user,
+    the hard limits are restrictions on the value imposed by the model.
+    For example, the thickness of a layer must be zero or more.
+
+    Any additional keyword arguments are preserved as properties of the
+    parameter. For example, *tip* and *units* for decorating an input form
+    in the GUI:
 
+         p = Parameter(10, name="width", units="cm", tip="Width of sample")
+
+    """
     # Parameters may be dependent on other parameters, and the
     # fit engine will need to access them.
+    #prior: Optional[BoundsType]
+    limits: Tuple[float, float]
+    _fixed: bool
+
     def parameters(self):
-        return [self]
+        pars = [self]
+        if hasattr(self.slot, 'parameters'):
+            pars += self.slot.parameters()
+        return pars
 
     def pmp(self, plus, minus=None, limits=None):
         """
@@ -100,7 +342,8 @@ def pmp(self, plus, minus=None, limits=None):
         The resulting range is converted to "nice" numbers.
         """
         bounds = mbounds.pmp(self.value, plus, minus, limits=limits)
-        self.bounds = mbounds.Bounded(*bounds)
+        self.bounds = bounds
+        self.fixed = False
         return self
 
     def pm(self, plus, minus=None, limits=None):
@@ -120,7 +363,8 @@ def pm(self, plus, minus=None, limits=None):
         The resulting range is converted to "nice" numbers.
         """
         bounds = mbounds.pm(self.value, plus, minus, limits=limits)
-        self.bounds = mbounds.Bounded(*bounds)
+        self.bounds = bounds
+        self.fixed = False
         return self
 
     def dev(self, std, mean=None, limits=None, sigma=None, mu=None):
@@ -145,25 +389,29 @@ def dev(self, std, mean=None, limits=None, sigma=None, mu=None):
                 mean = mu
         if mean is None:
             mean = self.value  # Note: value is an attribute of the derived class
-        if limits is None:
-            self.bounds = mbounds.Normal(mean, std)
-        else:
-            self.bounds = mbounds.BoundedNormal(mean, std, limits)
+        self.bounds = limits if limits is not None else (-inf, inf)
+        self.distribution = Normal(mean=mean, std=std)
+        self.fixed = False
         return self
 
-    def pdf(self, dist):
-        """
-        Allow the parameter to vary according to any continuous scipy.stats
-        distribution.
-        """
-        self.bounds = mbounds.Distribution(dist)
-        return self
+    # def pdf(self, dist):
+    #     """
+    #     Allow the parameter to vary according to any continuous scipy.stats
+    #     distribution.
+    #     """
+    #     # TODO: have to make some kind of registry for distributions?
+    #     # this will not work in new system of setting priors in model_reset.
+    #     self._set_bounds((-inf, inf))
+    #     self.distribution = dist
+    #     return self
 
     def range(self, low, high):
         """
         Allow the parameter to vary within the given range.
         """
-        self.bounds = mbounds.init_bounds((low, high))
+        self.bounds = (low, high)
+        self.distribution = Uniform()
+        self.fixed = False
         return self
 
     def soft_range(self, low, high, std):
@@ -171,96 +419,66 @@ def soft_range(self, low, high, std):
         Allow the parameter to vary within the given range, or with Gaussian
         probability, stray from the range.
         """
-        self.bounds = mbounds.SoftBounded(low, high, std)
+        self.bounds = (low, high)
+        self.distribution = UniformSoftBounded(std=std)
+        self.fixed = False
         return self
 
+    # Delegate to slots
+    @property
+    def value(self):
+        return int(self.slot) if self.discrete else float(self.slot)
+    @value.setter
+    def value(self, update):
+        self.slot.value = round(update) if self.discrete else update
     @property
-    def bounds(self):
-        """Fit bounds"""
-        # print "getting bounds for",self,self._bounds
-        return self._bounds
-
-    @bounds.setter
-    def bounds(self, b):
-        # print "setting bounds for",self
+    def fittable(self):
+        return isinstance(self.slot, Variable)
+    @property
+    def fixed(self):
+        return not self.fittable or self._fixed
+    @fixed.setter
+    def fixed(self, state):
+        # Can't set fixed to false if the parameter is not fittable
         if self.fittable:
-            self.fixed = (b is None)
-        self._bounds = b
+            self._fixed = state
+        elif not state:
+            raise TypeError(f"value in {self.name} is not fittable")
+
+    ## Use the following if bounds are on the value rather than the parameter
+    #@property
+    #def bounds(self):
+    #    return getattr(self.slot, 'bounds', None)
+    #@bounds.setter
+    #def bounds(self, b):
+    #    if not hasattr(self.slot, 'bounds'):
+    #        raise TypeError(f"{self.name} is not fittable so bounds can't be set")
+    #    if self.slot.fittable:
+    #        self.slot.fixed = (b is None)
+    #    self.slot.bounds = b
 
     # Functional form of parameter value access
     def __call__(self):
         return self.value
 
-    # Parameter algebra: express relationships between parameters
-    def __gt__(self, other):
-        return Constraint(self, other, "GT", ">")
-
-    def __ge__(self, other):
-        return Constraint(self, other, "GE", ">=")
-    def __le__(self, other):
-        return Constraint(self, other, "LE", "<=")
-
-    def __lt__(self, other):
-        return Constraint(self, other, "LT", "<")
-
-    # def __eq__(self, other):
-    #     return ConstraintEQ(self, other)
-
-    # def __ne__(self, other):
-    #     return ConstraintNE(self, other)
-
-    def __add__(self, other):
-        return Operator(self, other, "add", "+")
-
-    def __sub__(self, other):
-        return Operator(self, other, "sub", "-")
-
-    def __mul__(self, other):
-        return Operator(self, other, "mul", "*")
-
-    def __div__(self, other):
-        return Operator(self, other, "truediv", "/")
-
-    def __pow__(self, other):
-        return Operator(self, other, "pow", "**")
-
-    def __radd__(self, other):
-        return Operator(other, self, "add", "+")
-
-    def __rsub__(self, other):
-        return Operator(other, self, "sub", "-")
-
-    def __rmul__(self, other):
-        return Operator(other, self, "mul", "*")
-
-    def __rdiv__(self, other):
-        return Operator(other, self, "truediv", "/")
-
-    def __rpow__(self, other):
-        return Operator(other, self, "pow", "**")
-
-    def __abs__(self):
-        return _abs(self)
-
-    def __neg__(self):
-        return self * -1
-
-    def __pos__(self):
-        return self
-
     def __float__(self):
         return float(self.value)
 
-    __truediv__ = __div__
-    __rtruediv__ = __rdiv__
-
-    def nllf(self):
+    def nllf(self) -> float:
         """
         Return -log(P) for the current parameter value.
         """
-        return self.bounds.nllf(self.value)
+        value = self.value
+        if not (self.limits[0] <= value <= self.limits[1]):
+            # quick short-circuit if not meeting own limits:
+            return np.inf
+        else:
+            logp = self.prior.nllf(value)
+            if hasattr(self.slot, 'nllf'):
+                logp += self.slot.nllf()
+            return logp
 
-    def residual(self):
+    def residual(self) -> float:
         """
         Return the z score equivalent for the current parameter value.
 
@@ -272,7 +490,7 @@ def residual(self):
         compute the cdf of value in the parameter distribution and invert
         it using the ppf from the standard normal distribution.
         """
-        return self.bounds.residual(self.value)
+        return 0.0 if self.prior is None else self.prior.residual(self.value)
 
     def valid(self):
         """
@@ -284,7 +502,7 @@ def format(self):
         """
         Format the parameter, value and range as a string.
         """
-        return "%s=%g in %s" % (self, self.value, self.bounds)
+        return "%s=%g in %s" % (self, self.value, self.prior)
 
     def __str__(self):
         name = self.name if self.name is not None else '?'
@@ -293,83 +511,19 @@ def __str__(self):
     def __repr__(self):
         return "Parameter(%s)" % self
 
-    def to_dict(self):
-        """
-        Return a dict represention of the object.
-        """
-        # When reconstructing a model from json we will need to tie parameters
-        # together that were tied before. This can be done by managing a
-        # cache of allocated parameters indexed by id, and pulling from that
-        # cache on recontruction if the id already exists, otherwise create
-        # a new entry. Conveniently, this will handle free variable references
-        # in parameter sets as well. Note that the entire parameter description
-        # will be repeated each time it occurs, but there should be few
-        # enough of these that it isn't a problem.
-        # TODO: use id that is stable from session to session.
-        # TODO: have mechanism for clearing cache between save/load.
-        return dict(
-            type=type(self).__name__,
-            id=id(self), # Warning: this will be different every session
-            name=self.name,
-            value=self.value,
-            fixed=self.fixed,
-            fittable=self.fittable,
-            bounds=to_dict(self._bounds),
-            )
-
-
-class Constant(BaseParameter):
-    """
-    An unmodifiable value.
-    """
-    fittable = False
-    fixed = True
-
-    @property
-    def value(self):
-        return self._value
-
-    def __init__(self, value, name=None):
-        self._value = value
-        self.name = name
-
-    # to_dict() can inherit from BaseParameter
-
-
-class Parameter(BaseParameter):
-    """
-    A parameter is a symbolic value.
-
-    It can be fixed or it can vary within bounds.
-
-    p = Parameter(3).pmp(10)    # 3 +/- 10%
-    p = Parameter(3).pmp(-5,10) # 3 in [2.85,3.3] rounded to 2 digits
-    p = Parameter(3).pm(2)      # 3 +/- 2
-    p = Parameter(3).pm(-1,2)   # 3 in [2,5]
-    p = Parameter(3).range(0,5) # 3 in [0,5]
-
-    It has hard limits on the possible values, and a range that should live
-    within those hard limits.  The value should lie within the range for
-    it to be valid.  Some algorithms may drive the value outside the range
-    in order to satisfy soft It has a value which should lie within the range.
-
-    Other properties can decorate the parameter, such as tip for tool tip
-    and units for units.
-    """
-    fittable = True
-
+    # TODO: deprecate
     @classmethod
-    def default(cls, value, **kw):
+    def default(
+            cls: type,
+            value: Union[float, Tuple[float, float], ValueType],
+            **kw) -> "Parameter":
         """
-        Create a new parameter with the *value* and *kw* attributes, or return
-        the existing parameter if *value* is already a parameter.
-
-        The attributes are the same as those for Parameter, or whatever
-        subclass *cls* of Parameter is being created.
+        Create a new parameter with the *value* and *kw* attributes. If value
+        is already a parameter or expression, set it to that value.
         """
         # Need to constrain the parameter to fit within fixed limits and
         # to receive a name if a name has not already been provided.
-        if isinstance(value, BaseParameter):
+        if isinstance(value, Parameter):
             return value
         else:
             return cls(value, **kw)
@@ -378,23 +532,37 @@ def set(self, value):
         """
         Set a new value for the parameter, ignoring the bounds.
         """
-        self.value = value
+        self.slot.value = value
 
     def clip_set(self, value):
         """
         Set a new value for the parameter, clipping it to the bounds.
         """
-        low, high = self.bounds.limits
-        self.value = min(max(value, low), high)
-
-    def __init__(self, value=None, bounds=None, fixed=None, name=None, **kw):
-        # UI nicities:
-        # 1. check if we are started with value=range or bounds=range; if we
+        low, high = self.prior.limits
+        self.slot.value = builtins.min(builtins.max(value, low), high)
+
+    def __init__(
+            self,
+            value: Optional[Union[float, Tuple[float, float]]] = None,
+            slot: Optional[Union['Variable', ValueType]] = None,
+            #bounds: Optional[Union[BoundsType, Tuple[float, float]]]=None,
+            fixed: Optional[bool]=None,
+            name: Optional[str]=None,
+            id: Optional[str]=None,
+            limits: Optional[Tuple[Union[float, Literal[None, "-inf"]], Union[float, Literal[None, "inf"]]]]=None,
+            bounds: Optional[Tuple[Union[float, Literal["-inf"]], Union[float, Literal["inf"]]]]=None,
+            distribution: DistributionType = Uniform(),
+            discrete: bool = False,
+            **kw):
+        # Check if we are started with value=range or bounds=range; if we
         # are given bounds, then assume this is a fitted parameter, otherwise
         # the parameter defaults to fixed; if value is not set, use the
         # midpoint of the range.
         if bounds is None:
             try:
+                # Note: throws TypeError if not a sequence (which we want to
+                # fall through to the remainder of the function), or ValueError
+                # if the sequence is the wrong length (which we want to fail).
                 lo, hi = value
                 warnings.warn(DeprecationWarning("parameters can no longer be initialized with a fit range"))
                 bounds = lo, hi
@@ -403,86 +571,472 @@ def __init__(self, value=None, bounds=None, fixed=None, name=None, **kw):
                 pass
         if fixed is None:
             fixed = (bounds is None)
-        bounds = mbounds.init_bounds(bounds)
-        if value is None:
-            value = bounds.start_value()
-
-        # Store whatever values the user needs to associate with the parameter
-        # Models should set units and tool tips so the user interface has
-        # something to work with.
-        limits = kw.get('limits', (-inf, inf))
-        for k, v in kw.items():
-            setattr(self, k, v)
+        if slot is None:
+            if value is None:
+                value = float(bounds[0]) if bounds is not None else 0 # ? what else to do here?
+            if isinstance(value, (float, int)):
+                value = round(value) if discrete else value
+                slot = Variable(value)
+            elif isinstance(value, ValueProtocol):
+                slot = value
+            else:
+                raise TypeError("value %s: %s cannot be converted to Variable" % (str(name), str(value)))
+        assert isinstance(slot, (float, Variable, Expression, Parameter, Constant, Calculation))
 
-        # Initialize bounds, with limits clipped to the hard limits for the
-        # parameter
-        def clip(x, a, b):
-            return min(max(x, a), b)
+        self.slot = slot
+        self.name = name
+        self.id = id if id is not None else str(uuid.uuid4())
+        if limits is None:
+            limits = (-np.inf, np.inf)
+        self.limits  = (
+            (-np.inf if limits[0] is None else float(limits[0])),
+            (np.inf if limits[1] is None else float(limits[1])))
+        if bounds is not None:
+            bounds = (float(bounds[0]), float(bounds[1]))
         self.bounds = bounds
-        self.bounds.limits = (clip(self.bounds.limits[0], *limits),
-                              clip(self.bounds.limits[1], *limits))
-        self.value = value
+        self.distribution = distribution
+        # Note: fixed is True unless fixed=False or bounds=bounds were given
+        # as function arguments. Note that _set_bounds() will always set the
+        # fixed to False, so we need to reset it after calling _set_bounds().
         self.fixed = fixed
-        self.name = name
+        self.discrete = discrete
+
+        # Store whatever values the user needs to associate with the parameter.
+        # For example, models can set units and tool tips so the user interface
+        # has something to work with.
+        for k, v in kw.items():
+            setattr(self, k, v)
+        self.prior = None # to be filled by model_reset
 
     def randomize(self, rng=None):
         """
         Set a random value for the parameter.
         """
-        self.value = self.bounds.random(rng if rng is not None else mbounds.RNG)
+        self.value = self.prior.random(rng if rng is not None else mbounds.RNG)
 
     def feasible(self):
         """
         Value is within the limits defined by the model
         """
-        return self.bounds.limits[0] <= self.value <= self.bounds.limits[1]
+        return self.prior.limits[0] <= self.value <= self.prior.limits[1]
+
+    def equals(self, expression: ValueType):
+        """
+        Set a parameter equal to another parameter or expression.
+
+        If *expression=None* then free the parameter by giving it is own
+        slot with value equal to the present value of the expression, and
+        its bounds.
+        """
+        if isinstance(self.slot, Calculation):
+            raise TypeError("parameter is calculated by the model and cannot be changed")
+        elif expression is self:
+            # don't make a circular reference to self.
+            warnings.warn(f"{self} tried to make circular reference to self...")
+            pass
+        else:
+            self.slot = expression
+    
+    def unlink(self):
+        if isinstance(self.slot, Calculation):
+            raise TypeError("parameter is calculated by the model and cannot be changed")
+        # Replace the slot with a new variable initialized to the only variable value
+        self.slot = Variable(self.value)
+
+    def __copy__(self):
+        """ copy will only be called when a new instance is desired, with a different id """
+        obj = type(self).__new__(self.__class__)
+        obj.__dict__.update(self.__dict__)
+        obj.id = str(uuid.uuid4())
+        return obj
+
+
+@dataclass
+class Variable(ValueProtocol):
+    """
+    Saved state for a random variable in the model.
+    """
+    value: float
+    def __init__(self, value):
+        self.value = value
+    def parameters(self):
+        return []
 
-    # to_dict() can inherit from BaseParameter
 
+@schema(frozen=True, init=False)
+class Constant(ValueProtocol): # type: ignore
+    """
+    Saved state for an unmodifiable value.
 
-class Reference(Parameter):
+    A constant is like a fixed parameter. You can't change it's value, set
+    it equal to another parameter, or assign a prior distribution.
     """
-    Create an adaptor so that a model attribute can be treated as if it
-    were a parameter.  This allows only direct access, wherein the
-    storage for the parameter value is provided by the underlying model.
-
-    Indirect access, wherein the storage is provided by the parameter, cannot
-    be supported since the parameter has no way to detect that the model
-    is asking for the value of the attribute.  This means that model
-    attributes cannot be assigned to parameter expressions without some
-    trigger to update the values of the attributes in the model.
+
+    value: float
+    name: Optional[str]
+    id: str = field(metadata={"format": "uuid"})
+
+    fittable = False  # class property fixed across all objects
+    fixed = True # class property fixed across all objects
+
+    def __init__(self, value: float, name: Optional[str]=None, id: Optional[str]=None):
+        object.__setattr__(self, "value", value)
+        object.__setattr__(self, "name", name)
+        _id = id if id is not None else str(uuid.uuid4())
+        object.__setattr__(self, "id", _id)
+
+    def parameters(self):
+        return [self]
+
+    def __str__(self):
+        return self.name
+
+# ==== Arithmetic operators ===
+class Operators(str, Enum):
+    """ Operators that can be used to construct Expressions """
+
+    # operators including abs() are defined in _build_operator_mixin()
+    # functions are defined in numpy or in UserFunction (for min/max)
+
+    # unary operator
+    neg = "neg"
+    pos = "pos"
+    # binary operator
+    add = "add"
+    sub = "sub"
+    mul = "mul"
+    truediv = "div"
+    div = "intdiv"
+    pow = "pow"
+    # unary functional
+    #float = "float"  => float makes values concrete
+    #int = "int"  => values must be float; use floor, trunc, ceil, round
+    abs = "abs"
+
+    # unary functions
+    exp = "exp"
+    expm1 = "expm1"
+    log = "log"
+    log10 = "log10"
+    log1p = "log1p"
+    sqrt = "sqrt"
+    degrees = "degrees"
+    radians = "radians"
+    sin = "sin"
+    cos = "cos"
+    tan = "tan"
+    arcsin = "arcsin"
+    arccos = "arccos"
+    arctan = "arctan"
+    sinh = "sinh"
+    cosh = "cosh"
+    tanh = "tanh"
+    arcsinh = "arcsinh"
+    arccosh = "arccosh"
+    arctanh = "arctanh"
+    ceil = "ceil"
+    floor = "floor"
+    trunc = "trunc"
+    rint = "rint"
+    round = "round"  # round(a) => rint(a)
+    # binary functions
+    arctan2 = "arctan2"
+
+    # n-ary
+    min = "min"  # from builtins
+    max = "max"  # from builtins
+    # TODO: support sum(seq) and prod(seq) for tuple and list
+
+# Precedence for the python operators as given in manual. Numbers start
+# from one at the bottom of the table. The value itself is "highest" precedence
+# with a value of zero.
+# https://docs.python.org/3/reference/expressions.html#operator-precedence
+VALUE_PRECEDENCE = 0
+CALL_PRECEDENCE = 2
+OPERATOR_PRECEDENCE = {
+    "pow": 4,
+    "pos": 5,
+    "neg": 5,
+    "mul": 6,
+    "div": 6,
+    "truediv": 6,
+    "intdiv": 6,
+    "add": 7,
+    "sub": 7,
+    "gt": 12,
+    "lt": 12,
+    "ge": 12,
+    "le": 12,
+    "eq": 12,
+    "ne": 12,
+}
+OPERATOR_STRING = {
+    "pow": "**",
+    "pos": "+",
+    "neg": "-",
+    "mul": "*",
+    "div": "/",
+    "truediv": "/",
+    "intdiv": "//",
+    "add": "+",
+    "sub": "-",
+    "gt": ">",
+    "lt": "<",
+    "ge": ">=",
+    "le": "<=",
+    "eq": "==",
+    "ne": "!=",
+}
+
+
+def _lookup_operator(op_name):
+    if not hasattr(Operators, op_name) and op_name not in UserFunction.registry:
+        raise ValueError(f"function {op_name} is not available")
+    fn = None
+    # Check plugins first so we can override lookups in operator and numpy.
+    # This is needed for min/max.
+    if fn is None: # plugin functions from UserFunction registry
+        fn = UserFunction.registry.get(op_name, None)
+    if fn is None:
+        fn = getattr(operator, op_name, None) # operators from operators
+    if fn is None: # math functions from numpy
+        fn = getattr(np, op_name, None)
+    if fn is None:
+        raise RuntimeError(f"should not be here: {op_name} not found")
+    return fn
+
+def _precedence(obj: Any) -> int:
     """
+    Return operator precedence according to the python parsing hierarchy.
 
-    def __init__(self, obj, attr, **kw):
-        self.obj = obj
-        self.attr = attr
-        kw.setdefault('name', ".".join([obj.__class__.__name__, attr]))
-        Parameter.__init__(self, **kw)
+    Lower values are higher precedence. Values start at 0 for constants and
+    variables, and go up from there. Not all operators are covered.
+    """
+    if isinstance(obj, Expression):
+        return OPERATOR_PRECEDENCE.get(obj.op.name, CALL_PRECEDENCE)
+    return VALUE_PRECEDENCE
+
+@schema(init=False, eq=False, frozen=True)
+class Expression(ValueProtocol):
+    """
+    Parameter expression
+    """
+    fittable = False
+    fixed = True
+
+    op: Union[Operators, 'UserFunction'] # Enumerated str type {function_name: display_name}
+    args: Sequence[ValueType]
+    _fn: Callable[..., float] # _fn(float, float, ...) -> float
+
+    def __init__(self, op: Union[str, Operators, 'UserFunction'], args):
+        op = op if (isinstance(op, Operators) or isinstance(op, UserFunction)) else getattr(Operators, op)
+        object.__setattr__(self, 'op', op)
+        object.__setattr__(self, '_fn', _lookup_operator(op.name))
+        object.__setattr__(self, 'args', args)
+
+    def parameters(self):
+        # Walk expression tree combining parameters from each subexpression
+        return sum(
+            (v.parameters() for v in self.args if hasattr(v, 'parameters')),
+            [])
 
     @property
     def value(self):
-        return getattr(self.obj, self.attr)
+        return self._fn(*(float(arg) for arg in self.args))
 
-    @value.setter
-    def value(self, value):
-        setattr(self.obj, self.attr, value)
+    @property
+    def name(self):
+        return str(self)
 
-    def to_dict(self):
-        ret = Parameter.to_dict(self)
-        ret["attr"] = self.attr
-        # TODO: another impossibility---an arbitrary python object
-        # Clearly we need a (safe??) json pickler to handle the full
-        # complexity of an arbitrary model.
-        ret["obj"] = to_dict(self.obj)
-        return ret
+    def __str__(self):
+        prec = _precedence(self)
+        vals = [str(v) for v in self.args]
+        if self.op.name in ("pos", "neg"):
+            # +- a with parens as needed
+            a = f"({vals[0]})" if prec < _precedence(self.args[0]) else vals[0]
+            return f"{OPERATOR_STRING[self.op.name]}{a}"
+        elif self.op.name in ("add", "sub", "mul", "div", "truediv", "pow"):
+            # a +-*/** b with parens as needed
+            a = f"({vals[0]})" if prec < _precedence(self.args[0]) else vals[0]
+            b = f"({vals[1]})" if prec < _precedence(self.args[1]) else vals[1]
+            return f"{a} {OPERATOR_STRING[self.op.name]} {b}"
+        else:
+            # f(a, b, ...) with no parens needed
+            return f"{self.op.name}({', '.join(v for v in vals)})"
+
+def _make_unary_op(op_name: str):
+    op = getattr(Operators, op_name)
+    # Note: self is Parameter or Expression
+    fn = lambda self: Expression(op, (self,))
+    setattr(OperatorMixin, f'__{op_name}__', fn)
+
+def _make_binary_op(op_name: str):
+    op = getattr(Operators, op_name)
+    def fn(self, other):
+        return Expression(op, (self, other))
+    setattr(OperatorMixin, f'__{op_name}__', fn)
+    def rfn(self, other):
+        return Expression(op, (other, self))
+    setattr(OperatorMixin, f'__r{op_name}__', rfn)
+
+def _make_math_fn(fn_name: str):
+    op = getattr(Operators, fn_name)
+    def fn(*args): # first of args is self
+        if any([isinstance(arg, ValueProtocol) for arg in args]):
+            return Expression(op, args)
+        else:
+            # then all the args are floats: just return a float!
+            realized_fn = _lookup_operator(op.name)
+            return realized_fn(*args)
+    # define sin, etc., in the parameter and expression so that np.sin(a)
+    # will resolve to Expression('sin', tuple(a)), etc.
+    setattr(OperatorMixin, fn_name, fn)
+    # The np.sin(a) trick only works for a limited set of functions
+    # defined by numpy itself. For arbitrary user defined functions
+    # we add them to the bumps.pmath namespace so the user can find them.
+    setattr(pmath, fn_name, fn)
+
+def _build_operator_mixin():
+    unary_op = set(("pos", "neg", "abs"))
+    binary_op = set(("add", "sub", "mul", "div", "truediv", "pow"))
+    math_fn = set(v.name for v in Operators) - unary_op - binary_op
+    for op_name in unary_op:
+        _make_unary_op(op_name)
+    for op_name in binary_op:
+        _make_binary_op(op_name)
+    # By adding the math functions to the mixin, calling np.sin(parameter) or
+    # np.sin(expression) will return the generated expression for the object.
+    for fn_name in math_fn:
+        _make_math_fn(fn_name)
+_build_operator_mixin()
+
+# TODO: allow schema validation on user-defined functions
+class UserFunction:
+    """
+    User-defined functions.
+
+    This is a helper class for the @function decorator, which treats the
+    operator as one of the possible expression operators.
+
+    These won't be properly serialized/deserialized through the JSON schema
+    unless the function is registered in advance. The schema will not include
+    these functions as possible values even if registered, so a schema
+    validator may fail on one of these functions.
+    """
+    name: str
+    fn: Callable[..., float]
+    # A function registry to remember the code associated with the name.
+    # This is a class attribute, so it is initialized with an empty dict().
+    # Ignore complaints from lint.
+    # TODO: use pmath as our registry of available functions.
+    registry: Dict[str, Callable[..., float]] = {}
+    def __init__(self, fn: Callable):
+        name = fn.__name__
+        if name in UserFunction.registry:
+            raise TypeError(f"Function {name} already registered in bumps.")
+        UserFunction.registry[name] = fn
+        self.name = name
+
+def function(fn: Callable):
+    """
+    Convert a function into a delayed evaluator.
+
+    The value of the function is computed from the values of the parameters
+    at the time that the function value is requested rather than when the
+    function is created.
+    """
+    name = fn.__name__
+    op = UserFunction(fn)
+    def wrapped(*args: 'ValueType'):
+        return Expression(op, args)
+    wrapped.__name__ = fn.__name__
+    wrapped.__doc__ = fn.__doc__ if fn.__name__.endswith("d") else f"{fn.__name__}(Parameter)"
+    # Add the symbol to pmath
+    setattr(pmath, name, wrapped)
+    pmath.__all__.append(name)
+    return wrapped
+
+# min/max
+min = function(builtins.min)
+max = function(builtins.max)
+
+# Trig functions defined in degrees rather than radians.
+@function
+def cosd(v):
+    """Return the cosine of x (measured in in degrees)."""
+    return np.cos(np.radians(v))
+
+@function
+def sind(v):
+    """Return the sine of x (measured in in degrees)."""
+    return np.sin(np.radians(v))
+
+@function
+def tand(v):
+    """Return the tangent of x (measured in in degrees)."""
+    return np.tan(np.radians(v))
+
+@function
+def arccosd(v):
+    """Return the arc cosine (measured in in degrees) of x."""
+    return np.degrees(np.arccos(v))
 
+@function
+def arcsind(v):
+    """Return the arc sine (measured in in degrees) of x."""
+    return np.degrees(np.arcsin(v))
+
+@function
+def arctand(v):
+    """Return the arc tangent (measured in in degrees) of x."""
+    return np.degrees(np.arctan(v))
 
-class ParameterSet(object):
+@function
+def arctan2d(dy, dx):
+    """Return the arc tangent (measured in in degrees) of y/x.
+    Unlike atan(y/x), the signs of both x and y are considered."""
+    return np.degrees(np.arctan2(dy, dx))
+
+# Aliases for arcsin, etc., both here in bumps.parameters and in bumps.pmath.
+pmath.asin = asin = pmath.arcsin
+pmath.acos = acos = pmath.arccos
+pmath.atan = atan = pmath.arctan
+pmath.atan2 = atan2 = pmath.arctan2
+
+pmath.asind = asind = arcsind
+pmath.acosd = acosd = arccosd
+pmath.atand = atand = arctand
+pmath.atan2d = atan2d = arctan2d
+
+pmath.asinh = asinh = pmath.arcsinh
+pmath.acosh = acosh = pmath.arccosh
+pmath.atanh = atanh = pmath.arctanh
+
+pmath.__all__.extend((
+    'asin', 'acos', 'atan', 'atan2',
+    'asind', 'acosd', 'atand', 'atan2d',
+    'asinh', 'acosh', 'atanh',
+    ))
+
+#restate these for export, now that they're all defined:
+ValueType = Union[Parameter, Expression, Calculation, float]
+
+@schema(classname="ParameterSet")
+class ParameterSetSchema:
     """
     A parameter that depends on the model.
     """
+    names: Optional[List[str]]
+    reference: Parameter
+    parameterlist: Optional[List[Parameter]]
+
+class ParameterSet(ParameterSetSchema):
 
-    def __init__(self, reference, names=None):
+    def __init__(self,
+            reference: Parameter,
+            names: Optional[List[str]] = None,
+            parameterlist: Optional[List[Parameter]] = None
+        ):
         """
         Create a parameter set, with one parameter for each model name.
 
@@ -493,17 +1047,28 @@ def __init__(self, reference, names=None):
 
         *parameters* will be created, with one parameter per model.
         """
+        names = names if names is not None else []
         self.names = names
         self.reference = reference
         # TODO: explain better why parameters are using np.array
         # Force numpy semantics on slice operations by using an array
         # of objects rather than a list of objects
-        self.parameters = np.array([copy(reference) for _ in names])
+        if parameterlist is not None:
+            # we are being reinitialized with parameters
+            self.parameters = np.array(parameterlist)
+        else:
+            self.parameters = np.array([copy(reference) for _ in names])
         # print self.reference, self.parameters
         for p, n in zip(self.parameters, names):
             p.name = " ".join((n, p.name))
-        # Reference is no longer directly fittable
-        self.reference.fittable = False
+
+        # N.B. if the reference parameter is not referenced anywhere in the models,
+        # it will no longer show up in FitProblem.parameters
+        #self.__class__.parameterlist = property(self._get_parameterlist) #lambda self: self.parameters.tolist())
+
+    @property
+    def parameterlist(self) -> List[Parameter]:
+        return self.parameters.tolist()
 
     def to_dict(self):
         return {
@@ -598,6 +1163,7 @@ def pmp(self, *args, **kw):
             p.pmp(*args, **kw)
 
 
+@schema(init=False)
 class FreeVariables(object):
     """
     A collection of parameter sets for a group of models.
@@ -619,14 +1185,21 @@ class FreeVariables(object):
     new parameters for each model for the free parameters.  Setting up
     these copies was inconvenient.
     """
+    names: List[str]
+    parametersets: Dict[str, ParameterSet]
 
-    def __init__(self, names=None, **kw):
+    def __init__(self, names=None, parametersets=None, **kw):
         if names is None:
             raise TypeError("FreeVariables needs name=[model1, model2, ...]")
         self.names = names
-
-        # Create slots to hold the free variables
-        self._parametersets = dict((k, ParameterSet(v, names=names))
+        if parametersets is not None:
+            # assume that we are initializing with a dict of
+            # fully initialized ParameterSet objects
+            self.parametersets = parametersets
+        else:
+            # we are initializing with kw = Dict[key, (list of Parameters)]
+            # Create slots to hold the free variables
+            self.parametersets = dict((k, ParameterSet(v, names=names))
                                    for k, v in kw.items())
 
     # Shouldn't need explicit __getstate__/__setstate__ but mpi4py pickle
@@ -642,7 +1215,7 @@ def __getattr__(self, k):
         Return the parameter set for the given free parameter.
         """
         try:
-            return self._parametersets[k]
+            return self.parametersets[k]
         except KeyError:
             raise AttributeError('FreeVariables has no attribute %r' % k)
 
@@ -650,285 +1223,27 @@ def parameters(self):
         """
         Return the set of free variables for all the models.
         """
-        return dict((k, v.parameters) for k, v in self._parametersets.items())
+        return dict((k, v.parameters) for k, v in self.parametersets.items())
 
     def to_dict(self):
         return {
             'type': type(self).__name__,
             'names': self.names,
-            'parameters': to_dict(self._parametersets)
+            'parameters': to_dict(self.parametersets)
         }
 
     def set_model(self, i):
         """
         Set the reference parameters for model *i*.
         """
-        for p in self._parametersets.values():
+        for p in self.parametersets.values():
             p.set_model(i)
 
     def get_model(self, i):
         """
         Get the parameters for model *i* as {reference: substitution}
         """
-        return dict(p.get_model(i) for p in self._parametersets.values())
-
-# Current implementation computes values on the fly, so you only
-# need to plug the values into the parameters and the parameters
-# are automatically updated.
-#
-# This will not work well for wrapped models.  In those cases you
-# want to do a number of optimizations, such as only updating the
-#
-
-# ==== Comparison operators ===
-COMPARISONS = [
-    ('GT', '>'),
-    ('GE', '>='),
-    ('LE', '<='),
-    ('LT', '<'),
-    ('EQ', '=='),
-    ('NE', '!=')
-]
-
-class Constraint(object):
-    def __init__(self, a, b, op_name, op_str=""):
-        import operator
-        self.a, self.b = a, b
-        self.op_name = op_name
-        self.op = getattr(operator, op_name.lower())
-        self.op_str = op_str
-    def __bool__(self):
-        return self.op(float(self.a), float(self.b))
-    __nonzero__ = __bool__
-    def __str__(self):
-        return "(%s %s %s)" %(self.a, self.op_str, self.b)
-
-# ==== Arithmetic operators ===
-ALLOWED_OPERATORS = ["add","sub","mul","truediv","floordiv","pow"]
-
-class Operator(BaseParameter):
-    """
-    Parameter operator
-    """
-    def __init__(self, a, b, op_name, op_str):
-        import operator
-        if not op_name.lower() in ALLOWED_OPERATORS:
-            raise ValueError("Operator name %s is not in allowed operators: %s" % (op_name, str(ALLOWED_OPERATORS)))
-        self.a, self.b = a,b
-        self.op_name = op_name
-        self.op = getattr(operator, op_name.lower())
-        self.op_str = op_str
-        pars = []
-        if isinstance(a,BaseParameter): pars += a.parameters()
-        if isinstance(b,BaseParameter): pars += b.parameters()
-        self._parameters = pars
-        self.name = str(self)
-    def parameters(self):
-        return self._parameters
-    def to_dict(self):
-        return dict(
-            type="Operator",
-            op_name=self.op_name,
-            op_str=self.op_str,
-            left=to_dict(self.a),
-            right=to_dict(self.b),
-        )
-    @property
-    def value(self):
-        return self.op(float(self.a), float(self.b))
-    @property
-    def dvalue(self):
-        return float(self.a)
-    def __str__(self):
-        return "(%s %s %s)" % (self.a,self.op_str, self.b)
-
-def substitute(a):
-    """
-    Return structure a with values substituted for all parameters.
-
-    The function traverses lists, tuples and dicts recursively.  Things
-    which are not parameters are returned directly.
-    """
-    if isinstance(a, BaseParameter):
-        return float(a.value)
-    elif isinstance(a, tuple):
-        return tuple(substitute(v) for v in a)
-    elif isinstance(a, list):
-        return [substitute(v) for v in a]
-    elif isinstance(a, dict):
-        return dict((k, substitute(v)) for k, v in a.items())
-    elif isinstance(a, np.ndarray):
-        return np.array([substitute(v) for v in a])
-    else:
-        return a
-
-
-class Function(BaseParameter):
-    """
-    Delayed function evaluator.
-
-    f.value evaluates the function with the values of the
-    parameter arguments at the time f.value is referenced rather
-    than when the function was invoked.
-    """
-    __slots__ = ['op', 'args', 'kw']
-
-    def __init__(self, op, *args, **kw):
-        self.name = kw.pop('name', None)
-        self.op, self.args, self.kw = op, args, kw
-        self._parameters = self._find_parameters()
-
-    def _find_parameters(self):
-        # Figure out which arguments to the function are parameters
-        #deps = [p for p in self.args if isinstance(p,BaseParameter)]
-        args = [arg for arg in self.args if isinstance(arg, BaseParameter)]
-        kw = dict((name, arg) for name, arg in self.kw.items()
-                  if isinstance(arg, BaseParameter))
-        deps = flatten((args, kw))
-        # Find out which other parameters these parameters depend on.
-        res = []
-        for p in deps:
-            res.extend(p.parameters())
-        return res
-
-    def parameters(self):
-        return self._parameters
-
-    def _value(self):
-        # Expand args and kw, replacing instances of parameters
-        # with their values
-        return self.op(*substitute(self.args), **substitute(self.kw))
-    value = property(_value)
-
-    def to_dict(self):
-        return {
-            "type": "Function",
-            "name": self.name,
-            # TODO: function not stored properly in json
-            "op": to_dict(self.op),
-            "args": to_dict(self.args),
-            "kw": to_dict(self.kw),
-        }
-
-    def __getstate__(self):
-        return self.name, self.op, self.args, self.kw
-
-    def __setstate__(self, state):
-        self.name, self.op, self.args, self.kw = state
-        self._parameters = self._find_parameters()
-
-    def __str__(self):
-        if self.name is not None:
-            name = self.name
-        else:
-            args = [str(v) for v in self.args]
-            kw = [str(k) + "=" + str(v) for k, v in self.kw.items()]
-            name = self.op.__name__ + "(" + ", ".join(args + kw) + ")"
-        return name
-        #return "%s:%g" % (name, self.value)
-
-
-def function(op):
-    """
-    Convert a function into a delayed evaluator.
-
-    The value of the function is computed from the values of the parameters
-    at the time that the function value is requested rather than when the
-    function is created.
-    """
-    # Note: @functools.wraps(op) does not work with numpy ufuncs
-    # Note: @decorator does not work with builtins like abs
-    def function_generator(*args, **kw):
-        return Function(op, *args, **kw)
-    function_generator.__name__ = op.__name__
-    function_generator.__doc__ = op.__doc__
-    return function_generator
-_abs = function(abs)
-
-# Numpy trick: math functions from numpy delegate to the math function of
-# the class if that function exists as a class attribute.
-BaseParameter.exp = function(math.exp)
-BaseParameter.expm1 = function(math.expm1)
-BaseParameter.log = function(math.log)
-BaseParameter.log10 = function(math.log10)
-BaseParameter.log1p = function(math.log1p)
-BaseParameter.sqrt = function(math.sqrt)
-
-BaseParameter.degrees = function(math.degrees)
-BaseParameter.radians = function(math.radians)
-
-BaseParameter.sin = function(math.sin)
-BaseParameter.cos = function(math.cos)
-BaseParameter.tan = function(math.tan)
-BaseParameter.arcsin = function(math.asin)
-BaseParameter.arccos = function(math.acos)
-BaseParameter.arctan = function(math.atan)
-
-BaseParameter.sinh = function(math.sinh)
-BaseParameter.cosh = function(math.cosh)
-BaseParameter.tanh = function(math.tanh)
-BaseParameter.arcsinh = function(math.asinh)
-BaseParameter.arccosh = function(math.acosh)
-BaseParameter.arctanh = function(math.atanh)
-
-BaseParameter.ceil = function(math.ceil)
-BaseParameter.floor = function(math.floor)
-BaseParameter.trunc = function(math.trunc)
-
-def boxed_function(f):
-    box = function(f)
-    @wraps(f)
-    def wrapped(*args, **kw):
-        if any(isinstance(v, BaseParameter) for v in args):
-            return box(*args, **kw)
-        else:
-            return f(*args, **kw)
-    return wrapped
-
-# arctan2 is special since either argument can be a parameter
-arctan2 = boxed_function(math.atan2)
-
-# Trig functions defined in degrees rather than radians
-@boxed_function
-def cosd(v):
-    """Return the cosine of x (measured in in degrees)."""
-    return math.cos(math.radians(v))
-
-@boxed_function
-def sind(v):
-    """Return the sine of x (measured in in degrees)."""
-    return math.sin(math.radians(v))
-
-@boxed_function
-def tand(v):
-    """Return the tangent of x (measured in in degrees)."""
-    return math.tan(math.radians(v))
-
-@boxed_function
-def acosd(v):
-    """Return the arc cosine (measured in in degrees) of x."""
-    return math.degrees(math.acos(v))
-arccosd = acosd
-
-@boxed_function
-def asind(v):
-    """Return the arc sine (measured in in degrees) of x."""
-    return math.degrees(math.asin(v))
-arcsind = asind
-
-@boxed_function
-def atand(v):
-    """Return the arc tangent (measured in in degrees) of x."""
-    return math.degrees(math.atan(v))
-arctand = atand
-
-@boxed_function
-def atan2d(dy, dx):
-    """Return the arc tangent (measured in in degrees) of y/x.
-    Unlike atan(y/x), the signs of both x and y are considered."""
-    return math.degrees(math.atan2(dy, dx))
-arctan2d = atan2d
-
+        return dict(p.get_model(i) for p in self.parametersets.values())
 
 
 def flatten(s):
@@ -938,7 +1253,7 @@ def flatten(s):
         raise TypeError("parameter flattening cannot order sets")
     elif isinstance(s, dict):
         return reduce(lambda a, b: a + flatten(s[b]), sorted(s.keys()), [])
-    elif isinstance(s, BaseParameter):
+    elif isinstance(s, ValueProtocol):
         return [s]
     elif s is None:
         return []
@@ -946,13 +1261,14 @@ def flatten(s):
         raise TypeError("don't understand type %s for %r" % (type(s), s))
 
 
-def format(p, indent=0, freevars={}, field=None):
+def format(p, indent=0, freevars=None, field=None):
     """
     Format parameter set for printing.
 
     Note that this only says how the parameters are arranged, not how they
     relate to each other.
     """
+    freevars = {} if freevars is None else freevars
     p = freevars.get(id(p), p)
     if isinstance(p, dict) and p != {}:
         res = []
@@ -986,14 +1302,14 @@ def format(p, indent=0, freevars={}, field=None):
             s += str(p) + " = "
         s += "%g" % p.value
         if not p.fixed:
-            s += " in [%g,%g]" %  tuple(p.bounds.limits)
+            s += " in [%g,%g]" %  tuple(p.bounds)
         return s
 
-    elif isinstance(p, BaseParameter):
+    elif isinstance(p, Parameter):
         return "%s = %g" % (str(p), p.value)
 
     else:
-        return "None"
+        return str(p)
 
 
 def summarize(pars, sorted=False):
@@ -1005,12 +1321,12 @@ def summarize(pars, sorted=False):
     """
     output = []
     if sorted:
-        pars = sorted(pars, cmp=lambda x, y: cmp(x.name, y.name))
+        pars = sorted(pars, key=lambda x: x.name)
     for p in pars:
         if not isfinite(p.value):
             bar = ["*invalid* "]
         else:
-            position = int(p.bounds.get01(p.value) * 9.999999999)
+            position = int(p.prior.get01(p.value) * 9.999999999)
             bar = ['.'] * 10
             if position < 0:
                 bar[0] = '<'
@@ -1023,7 +1339,7 @@ def summarize(pars, sorted=False):
     return "\n".join(output)
 
 
-def unique(s):
+def unique(s) -> List[Parameter]:
     """
     Return the unique set of parameters
 
@@ -1058,10 +1374,10 @@ def fittable(s):
 
     Note that some fittable parameters may be fixed during the fit.
     """
-    return [p for p in unique(s) if not p.fittable]
+    return [p for p in unique(s) if p.fittable]
 
 
-def varying(s):
+def varying(s: List[Parameter]) -> List[Parameter]:
     """
     Return the list of fitted parameters in the model.
 
@@ -1069,32 +1385,113 @@ def varying(s):
     """
     return [p for p in unique(s) if not p.fixed]
 
+def _has_prior(p: Parameter) -> bool:
+    prior = getattr(p, 'prior', None)
+    limits = getattr(prior, 'limits', (-np.inf, np.inf))
+    return (
+        prior is not None
+        and not isinstance(prior, mbounds.Unbounded)
+        and limits != (-np.inf, np.inf))
+
+def priors(s: List[Parameter]) -> List[Parameter]:
+    """
+    Return the list of parameters (fitted or computed) that have prior
+    probabilities associated with them. This includes all varying parameters,
+    plus expressions (including simple links), but ignoring constants and
+    fixed parameters whose probabilities won't change the fits.
+    """
+    return [p for p in unique(s) if _has_prior(p)]
 
-def randomize(s):
+def randomize(s: List[Parameter]):
     """
     Set random values to the parameters in the parameter set, with
     values chosen according to the bounds.
     """
     for p in s:
-        p.value = p.bounds.random(1)[0]
+        p.value = p.prior.random(1)[0]
 
 
-def current(s):
+def current(s: List[Parameter]):
     return [p.value for p in s]
 
 # ========= trash ===================
 
+def copy_linked(has_parameters, free_names=None):
+    """
+    make a copy of an object with parameters
+     - then link all the parameters, except
+     - those with names matching "free_names"
+    """
+    assert callable(getattr(has_parameters, 'parameters', None)) == True
+    from copy import deepcopy
+    copied = deepcopy(has_parameters)
+    free_names = [] if free_names is None else free_names
+    original_pars = unique(has_parameters.parameters())
+    copied_pars = unique(copied.parameters())
+    for op, cp in zip(original_pars, copied_pars):
+        if not op.name in free_names:
+            cp.slot = op.slot
+        else:
+            cp.id = str(uuid.uuid4())
+    return copied
+
+# ==== Comparison operators ===
+class Comparisons(Enum):
+    """comparison operators"""
+
+    gt = '>'
+    ge = '>='
+    le = '<='
+    lt = '<'
+    #eq = '=='
+    #ne = '!='
+
 
-class IntegerParameter(Parameter):
-    discrete = True
+@schema(frozen=True)
+class Constraint:
+    """ Express inequality constraints between model elements """
+
+    fixed = True
+    
+    op: Comparisons
+    a: ValueType
+    b: ValueType
+
+    def __init__(self, a, b, op):
+        import operator
+        object.__setattr__(self, 'a', a)
+        object.__setattr__(self, 'b', b)
+        op_name = str(Comparisons(op).name)
+        object.__setattr__(self, 'compare', getattr(operator, op_name.lower()))
+        object.__setattr__(self, 'op', op)
+
+    # TODO: is this really necessary?  What is the reason for this trap?
+    # It seems like being able to cast with bool(Constraint) would be
+    # useful in some circumstances, like doing max(List[Parameter]), which
+    # currently fails.
+    def __bool__(self):
+        raise TypeError("failed bool")
+    __nonzero__ = __bool__
+    def __float__(self):
+        """return a float value that can be differentiated"""
+        return 0. if self.satisfied else abs(float(self.a) - float(self.b))
+    def __str__(self):
+        return "(%s %s %s)" %(self.a, self.op, self.b)
+    @property
+    def satisfied(self):
+        return self.compare(float(self.a), float(self.b))
 
-    def _get_value(self):
-        return self._value
 
-    def _set_value(self, value):
-        self._value = int(value)
-    value = property(_get_value, _set_value)
+def _make_constraint(op_str: str) -> Callable[..., Constraint]:
+    return lambda self, other: Constraint(self, other, op_str)
 
+def _build_constraints_mixin():
+    for comp_item in Comparisons:
+        op_name = comp_item.name
+        op_str = comp_item.value
+        setattr(OperatorMixin, f'__{op_name}__', _make_constraint(op_str))
+
+_build_constraints_mixin()
 
 class Alias(object):
     """
@@ -1105,7 +1502,6 @@ class Alias(object):
     parameter will have the full parameter semantics, including
     the ability to replace a fixed value with a parameter expression.
 
-    **Deprecated** :class:`Reference` does this better.
     """
 
     def __init__(self, obj, attr, p=None, name=None):
@@ -1130,10 +1526,233 @@ def to_dict(self):
             'attr': self.attr,
         }
 
+def substitute(a):
+    """
+    Return structure a with values substituted for all parameters.
+
+    The function traverses lists, tuples and dicts recursively.  Things
+    which are not parameters are returned directly.
+    """
+    if isinstance(a, ValueProtocol):
+        return float(a.value)
+    elif isinstance(a, tuple):
+        return tuple(substitute(v) for v in a)
+    elif isinstance(a, list):
+        return [substitute(v) for v in a]
+    elif isinstance(a, dict):
+        return dict((k, substitute(v)) for k, v in a.items())
+    elif isinstance(a, np.ndarray):
+        return np.array([substitute(v) for v in a])
+    else:
+        return a
+
+
+class Function(ValueProtocol):
+    """
+    **DEPRECATED**
+
+    Delayed function evaluator.
+
+    f.value evaluates the function with the values of the
+    parameter arguments at the time f.value is referenced rather
+    than when the function was invoked.
+    """
+    __slots__ = ['op', 'args', 'kw']
+    op: Callable[..., float]
+    args: Optional[Any]
+    kw: Dict[Any, Any]
+
+    def __init__(self, op, *args, **kw):
+        warnings.warn("Function no longer supported", DeprecationWarning, stacklevel=1)
+        self.name = kw.pop('name', None)
+        self.op, self.args, self.kw = op, args, kw
+        self._parameters = self._find_parameters()
+
+    def _find_parameters(self):
+        # Figure out which arguments to the function are parameters
+        #deps = [p for p in self.args if isinstance(p, ValueProtocol)]
+        args = [arg for arg in self.args if isinstance(arg, ValueProtocol)]
+        kw = dict((name, arg) for name, arg in self.kw.items()
+                  if isinstance(arg, ValueProtocol))
+        deps = flatten((args, kw))
+        # Find out which other parameters these parameters depend on.
+        res = []
+        for p in deps:
+            res.extend(p.parameters())
+        return res
+
+    def parameters(self):
+        return self._parameters
+
+    def _value(self):
+        # Expand args and kw, replacing instances of parameters
+        # with their values
+        return self.op(*substitute(self.args), **substitute(self.kw))
+    value = property(_value)
+
+    def to_dict(self):
+        return {
+            "type": "Function",
+            "name": self.name,
+            # TODO: function not stored properly in json
+            "op": to_dict(self.op),
+            "args": to_dict(self.args),
+            "kw": to_dict(self.kw),
+        }
+
+    def __getstate__(self):
+        return self.name, self.op, self.args, self.kw
+
+    def __setstate__(self, state):
+        self.name, self.op, self.args, self.kw = state
+        self._parameters = self._find_parameters()
+
+    def __str__(self):
+        if self.name is not None:
+            name = self.name
+        else:
+            args = [str(v) for v in self.args]
+            kw = [str(k) + "=" + str(v) for k, v in self.kw.items()]
+            name = self.op.__name__ + "(" + ", ".join(args + kw) + ")"
+        return name
+        #return "%s:%g" % (name, self.value)
+
+
+# ===== Tests ====
+
 def test_operator():
     a = Parameter(1, name='a')
     b = Parameter(2, name='b')
-    a_b = a + b
+    c = Parameter(3, name='c')
+    C = Constant(5, name='C')
+
+    assert a.fixed
+
+    # Check strings
+    assert str(a + b) == 'a + b'
+    assert (a+b).name == 'a + b'
+    assert str(-a) == '-a'
+    assert (-a).value == -a.value
+    assert str(a+b*c) == 'a + b * c'
+    assert str((a+b)*c) == '(a + b) * c'
+    assert str(np.sin(a+b)*c) == 'sin(a + b) * c'
+    assert str(a+C) == 'a + C'
+    assert str(a+C+3) == 'a + C + 3'
+    assert str(3 + a+C) == '3 + a + C'
+    assert str(a.sin()) == 'sin(a)'
+    assert str(atan2(a, b)) == 'arctan2(a, b)'
+    # float(expr) evaluates the expression; it doesn't build an expr with float.
+
+    # Check parameters
+    assert (a+b).parameters() == [a, b]
+    assert (np.sin(a+b)*c).parameters() == [a, b, c]
+
+    # Check values
     a.value = 3
-    assert a_b.value == 5.
-    assert a_b.name == '(a + b)'
\ No newline at end of file
+    assert (a+b).value == 5.
+    assert float(a+b) == a.value + b.value
+    assert a.sin().value == np.sin(a.value)
+    assert (3 + a+C).value == 3+3+5
+    assert np.sin(a+b).value == np.sin(a.value+b.value)
+    assert atan2(a, b).value == atan2(a.value, b.value)
+
+    # Make sure that evaluation is lazy. Capture the expression with one
+    # set of values for the parameters, update them with a new set of values,
+    # then check if the result is what you get when you call the function
+    # directly on those new values.
+    scope = locals() # record the currently available parameter handles
+    def capture_test(expr, result, **kw):
+        #print("checking", expr, "for", kw, "yields", result)
+        saved = {k: scope[k].value for k in kw}
+        for k, v in kw.items():
+            scope[k].value = float(v)
+        try:
+            assert expr.value == result, f"for {expr} expected {result} but got {expr.value}"
+        finally:
+            for k, v in saved.items():
+                scope[k].value = v
+    capture_test(np.sin(a+b), np.sin(0.5 + 3), a=0.5, b=3)
+    capture_test(np.arctan2(a, b), atan2(0.5, 3), a=0.5, b=3)
+    capture_test(np.round(a), np.round(-0.6), a=-0.6)
+    capture_test(min(a, b), builtins.min(-0.6, 3), a=-0.6, b=3)
+    capture_test(min(a, b, -2), builtins.min(-0.6, 3, -2), a=-0.6, b=3)
+    capture_test(abs(a), 2.5, a=-2.5)
+
+    # Check that symbols are defined in pmath
+    capture_test(pmath.sind(a), np.sin(np.radians(25)), a=25)
+    assert 'sind' in pmath.__all__
+
+    # TODO: can we evaluate an expression for an entire population at once?
+
+    # Check slots
+    limited = Parameter(3, name='limited', limits=[0.5, 1.5], bounds=[0, 1])
+    limited.add_prior()
+    assert np.isinf(limited.nllf())
+    assert np.isinf(limited.nllf())
+    limited.value = 0.6
+    assert limited.nllf() == 0.0
+    limited.value = 0.2
+    assert np.isinf(limited.nllf())
+
+    limited.equals(a+b)
+    assert limited.value == (a+b).value
+
+    assert np.isinf(limited.nllf())
+    a.value = b.value = 0.1
+    assert np.isinf(limited.nllf())
+    a.value = b.value = 0.3
+    assert limited.nllf() == 0.0
+    try:
+        limited.value = 5
+        failed = True
+    except Exception:
+        # TODO: define which error improper assignment should raise
+        # Currently this raises an attribute error on limited.slot.value
+        failed = False
+    if failed:
+        raise RuntimeError("failed to raise error when assigning value to expression")
+
+    # Check parameter list operations
+    s = [a, limited]
+    assert unique(s) == [a, limited, b]
+    assert fittable(s) == [a, b]
+    assert varying(s) == []
+    b.range(0, 3)
+    assert not b.fixed
+    assert varying(s) == [b]
+    assert current(s) == [a.value, limited.value]
+
+    # Check normal deviation
+    mu, sigma = 3, 2
+    b.dev(sigma, mean=mu)
+    b.value = 4
+    b.add_prior()
+    nllf_target = 0.5*((b.value-mu)/sigma)**2 + np.log(2*np.pi*sigma**2)/2
+    assert abs(b.nllf() - nllf_target)/nllf_target < 1e-12
+
+ # Check conditions
+    a.value, b.value = 3, 4
+    capture = a < b
+    assert isinstance(capture, Constraint)
+    assert capture.satisfied
+    a.value, b.value = 4, 3
+    assert not capture.satisfied
+
+    scope = locals()
+    def raises(condition_str, exception):
+        try:
+            eval(condition_str, locals=scope)
+        except exception:
+            pass
+        else:
+            raise AssertionError(f"{condition_str} does not raise {exception}")
+    raises("a < b < c", TypeError)
+    raises("a < b and b < c", TypeError)
+    raises("a < b or b < c", TypeError)
+    raises("not (a < b)", TypeError)
+    raises("not a", TypeError)
+    raises("a and b", TypeError)
+    raises("a or b", TypeError)
+
+if __name__ == "__main__":
+    test_operator()
diff --git a/bumps/plotutil.py b/bumps/plotutil.py
index 3f47e61e..e719e64f 100644
--- a/bumps/plotutil.py
+++ b/bumps/plotutil.py
@@ -38,7 +38,7 @@ def auto_shift(offset):
 
 # ======== Color functions ========
 
-def next_color():
+def next_color(axes=None):
     """
     Return the next color in the plot color cycle.
 
@@ -52,7 +52,9 @@ def next_color():
         plt.plot(x, y, '-', color=dhsv(color, dv=-0.2))
     """
     import pylab
-    lines = pylab.gca()._get_lines
+    if axes is None:
+        axes = pylab.gca()
+    lines = axes._get_lines
     try:
         base = next(lines.prop_cycler)['color']
     except Exception:
@@ -115,7 +117,7 @@ def dhsv(color, dh=0., ds=0., dv=0., da=0.):
 
 # ==== Specialized plotters =====
 
-def plot_quantiles(x, y, contours, color, alpha=None):
+def plot_quantiles(x, y, contours, color, alpha=None, axes=None):
     """
     Plot quantile curves for a set of lines.
 
@@ -134,18 +136,20 @@ def plot_quantiles(x, y, contours, color, alpha=None):
     default value, alpha=2./(#contours+1), works pretty well.
     """
     _, q = form_quantiles(y, contours)
-    _plot_quantiles(x, q, color, alpha)
+    _plot_quantiles(x, q, color, alpha, axes=axes)
 
-def _plot_quantiles(x, q, color, alpha):
+def _plot_quantiles(x, q, color, alpha, axes=None):
     import pylab
     # print "p",p
     # print "q",q[:,:,0]
     # print "y",y[:,0]
+    if axes is None:
+        axes = pylab.gca()
     if alpha is None:
         alpha = 2. / (len(q) + 1)
     edgecolor = dhsv(color, ds=-(1 - alpha), dv=(1 - alpha))
     for lo, hi in q:
-        pylab.fill_between(x, lo, hi,
+        axes.fill_between(x, lo, hi,
                            facecolor=color, edgecolor=edgecolor,
                            alpha=alpha)
 
diff --git a/bumps/plugin.py b/bumps/plugin.py
index 844b1cc8..e1e94ec0 100644
--- a/bumps/plugin.py
+++ b/bumps/plugin.py
@@ -83,7 +83,7 @@ def calc_errors(problem, sample):
     return None
 
 
-def show_errors(errs):
+def show_errors(errs, fig=None):
     """
     Display the model with uncertainty on the current figure.
 
diff --git a/bumps/pmath.py b/bumps/pmath.py
index e10a286a..131055e9 100644
--- a/bumps/pmath.py
+++ b/bumps/pmath.py
@@ -1,98 +1,23 @@
 """
 Standard math functions for parameter expressions.
 """
-import math
-
-from six.moves import reduce, builtins
-
-from .parameter import function
-__all__ = [
-    'exp', 'log', 'log10', 'sqrt',
-    'sin', 'cos', 'tan', 'asin', 'acos', 'atan', 'atan2',
-    'sind', 'cosd', 'tand', 'asind', 'acosd', 'atand', 'atan2d',
-    'sinh', 'cosh', 'tanh', 'asinh', 'acosh', 'atanh',
-    'degrees', 'radians',
-    'sum', 'prod',
-]
-
-def _cosd(v):
-    """Return the cosine of x (measured in in degrees)."""
-    return math.cos(math.radians(v))
-
-
-def _sind(v):
-    """Return the sine of x (measured in in degrees)."""
-    return math.sin(math.radians(v))
-
-
-def _tand(v):
-    """Return the tangent of x (measured in in degrees)."""
-    return math.tan(math.radians(v))
-
-
-def _acosd(v):
-    """Return the arc cosine (measured in in degrees) of x."""
-    return math.degrees(math.acos(v))
-
-
-def _asind(v):
-    """Return the arc sine (measured in in degrees) of x."""
-    return math.degrees(math.asin(v))
-
-
-def _atand(v):
-    """Return the arc tangent (measured in in degrees) of x."""
-    return math.degrees(math.atan(v))
-
-
-def _atan2d(dy, dx):
-    """Return the arc tangent (measured in in degrees) of y/x.
-    Unlike atan(y/x), the signs of both x and y are considered."""
-    return math.degrees(math.atan2(dy, dx))
-
-def _prod(s):
-    """Return the product of a sequence of numbers."""
-    return reduce(lambda x, y: x * y, s, 1)
-
-exp = function(math.exp)
-log = function(math.log)
-log10 = function(math.log10)
-sqrt = function(math.sqrt)
-
-degrees = function(math.degrees)
-radians = function(math.radians)
-
-sin = function(math.sin)
-cos = function(math.cos)
-tan = function(math.tan)
-asin = function(math.asin)
-acos = function(math.acos)
-atan = function(math.atan)
-atan2 = function(math.atan2)
-
-sind = function(_sind)
-cosd = function(_cosd)
-tand = function(_tand)
-asind = function(_asind)
-acosd = function(_acosd)
-atand = function(_atand)
-atan2d = function(_atan2d)
-
-sinh = function(math.sinh)
-cosh = function(math.cosh)
-tanh = function(math.tanh)
-asinh = function(math.asinh)
-acosh = function(math.acosh)
-atanh = function(math.atanh)
-
-sum = function(builtins.sum)
-prod = function(_prod)
-
-min = function(builtins.min)
-max = function(builtins.max)
-
-# Define pickler for numpy ufuncs
-#import copy_reg
-#def udump(f): return f.__name__
-#def uload(name): return getattr(np, name)
-#copy_reg.pickle(np.ufunc, udump, uload)
+__all__ = []
+
+# Note: the symbols in this module are defined dynamically by parameter.py
+
+# TODO: need a pmath test
+# Something like:
+#
+# import numpy as np
+# from bumps.parameter import Parameter, function
+#
+# # Define a plugin function
+# @function
+# def square(x): return x*x
+#
+# # Check all pmath symbols are imported, including plugin functions. This
+# # must be done after the plugin functions are registered.
+# from bumps.pmath import *
+# a = Parameter(value=0.2)
+# assert asind(a).value == np.degrees(np.arcsin(a.value))
+# assert square(a).value = a.value**2
diff --git a/bumps/serialize.py b/bumps/serialize.py
new file mode 100644
index 00000000..b6dcf799
--- /dev/null
+++ b/bumps/serialize.py
@@ -0,0 +1,194 @@
+import json
+
+from . import bounds, parameter, plugin
+
+import numpy as np
+from dataclasses import is_dataclass, fields, dataclass
+from importlib import import_module
+from typing import Dict, List, Literal, Optional, Tuple, Union
+from types import GeneratorType
+import traceback
+import warnings
+import asyncio
+from collections import defaultdict
+
+from .util import SCHEMA_ATTRIBUTE_NAME, schema, NumpyArray
+
+DEBUG = False
+MISSING = object()
+REFERENCE_TYPE = "Reference"
+
+@dataclass
+class Reference:
+    id: str
+    type: Literal[REFERENCE_TYPE]
+
+class Deserializer(object):
+    def __init__(self, loop=None):
+        self.loop = loop if loop is not None else asyncio.get_running_loop()
+        self.lock = asyncio.Lock()
+        self.deferred: Dict[str, asyncio.Future] = defaultdict(self.loop.create_future)
+    
+    def get_class(self, module_name: str, class_name: str):
+        return getattr(import_module(module_name), class_name)
+
+    async def rehydrate(self, obj):
+        if isinstance(obj, dict):
+            obj = obj.copy()
+            t: str = obj.pop('type', MISSING)
+            if t == REFERENCE_TYPE:
+                obj_id: str = obj.get("id", MISSING)
+                if obj_id is MISSING:
+                    raise ValueError("object id is required for Reference type")
+                async with self.lock:
+                    if DEBUG and not obj_id in self.deferred:
+                        print(f"creating deferred value: {obj_id}")
+                    fut = self.deferred[obj_id]
+                return await fut
+
+            if t == 'bumps.util.NumpyArray':
+                # skip processing values list in ndarray
+                return self.ndarray(obj)
+            else:
+                keys = obj.keys()
+                values = await asyncio.gather(*[self.rehydrate(value) for value in obj.values()])
+                for key,value in zip(keys, values):
+                    obj[key] = value
+
+            if t is MISSING:
+                # no "type" provided, so no class to instantiate
+                return obj
+            elif t == 'bumps.util.NumpyArray':
+                return self.ndarray(obj)
+            else:
+                try:
+                    module_name, class_name = t.rsplit('.', 1)
+                    klass = self.get_class(module_name, class_name)
+                    hydrated = await self.instantiate(klass, t, obj)
+                    return hydrated
+                except Exception:
+                    # there is a type, but it is not found...
+                    raise ValueError("type %s not found!" % t, obj)
+
+        elif isinstance(obj, list):
+            return await asyncio.gather(*[self.rehydrate(v) for v in obj])
+        elif isinstance(obj, tuple):
+            values = await asyncio.gather(*[self.rehydrate(v) for v in obj])
+            return tuple(values)
+        else:
+            return obj
+
+    async def instantiate(self, klass: type, typename: str, serialized: dict):
+        s = serialized.copy()
+        obj_id: str = s.get("id", MISSING) # will id be a required schema element?
+        #print('rehydrating: ', typename)
+        fut: Optional[asyncio.Future] = None
+        if obj_id is not MISSING:
+            async with self.lock:
+                # self.deferred is a defaultdict, so a new Future is returned if none
+                # already exists for that obj_id
+                fut = self.deferred[obj_id]
+                if fut.done():
+                    return fut.result()
+
+        # if klass provides 'from_dict' method, use it - 
+        # otherwise use klass.__init__ directly.
+        class_factory = getattr(klass, 'from_dict', klass)
+        try:
+            hydrated = class_factory(**s)
+        except Exception as e:
+            print(class_factory, s, typename)
+            raise e
+
+        if fut is not None:
+            if DEBUG:
+                print(f"setting future value for {obj_id} to {hydrated}")
+            async with self.lock:
+                fut.set_result(hydrated)
+            
+        return hydrated
+
+    def ndarray(self, obj: dict):
+        return np.asarray(obj['values'], dtype=np.dtype(obj.get('dtype', float)))
+
+
+async def async_from_dict(serialized, loop=None):
+    oasis = Deserializer(loop)
+    return await oasis.rehydrate(serialized)
+
+def from_dict(serialized):
+    return asyncio.run(async_from_dict(serialized))
+
+async def async_load(filename):
+    """ use in e.g. Jupyter notebooks """
+    with open(filename, 'r') as fid:
+        return await async_from_dict(json.loads(fid.read()))
+    
+def load(filename):
+    return asyncio.run(async_load(filename))
+
+class Serializer:
+    def __init__(self, use_refs=True):
+        self.refs = {}
+        self.use_refs = use_refs
+
+    def dataclass_to_dict(self, dclass, include=None, exclude=None):
+        all_fields = fields(dclass)
+        if include is not None:
+            all_fields = [f for f in all_fields if f.name in include]
+        elif exclude is not None:
+            all_fields = [f for f in all_fields if not f.name in exclude and not f.name.startswith('_')]
+        else:
+            all_fields = [f for f in all_fields if not f.name.startswith('_')]
+        cls = dclass.__class__
+        fqn = f"{cls.__module__}.{cls.__qualname__}"
+        output = dict([(f.name, self.to_dict(getattr(dclass, f.name))) for f in all_fields])
+        output["type"] = fqn
+        return output
+
+    def to_dict(self, obj):
+        if hasattr(obj, SCHEMA_ATTRIBUTE_NAME):
+            schema_opts = getattr(obj, SCHEMA_ATTRIBUTE_NAME)
+            include = schema_opts.get("include")
+            exclude = schema_opts.get("exclude")
+            if self.use_refs and hasattr(obj, 'id') and obj.id in self.refs:
+                return dict(id=obj.id, type=REFERENCE_TYPE)
+            else:
+                output = self.dataclass_to_dict(obj, include=include, exclude=exclude)
+                if self.use_refs and hasattr(obj, 'id'):
+                    self.refs.setdefault(obj.id, output)
+                return output
+        elif is_dataclass(obj):
+            return self.dataclass_to_dict(obj)
+        elif isinstance(obj, (list, tuple)):
+            return type(obj)(self.to_dict(v) for v in obj)
+        elif isinstance(obj, GeneratorType):
+            return list(self.to_dict(v) for v in obj)
+        # elif isinstance(obj, GeneratorType) and issubclass(obj_type, Tuple):
+        #     return tuple(to_dict(v) for v in obj)
+        elif isinstance(obj, dict):
+            return type(obj)((self.to_dict(k), self.to_dict(v))
+                            for k, v in obj.items())
+        elif isinstance(obj, np.ndarray) and obj.dtype.kind in ['f', 'i', 'U']:
+            return dict(type="bumps.util.NumpyArray", dtype=str(obj.dtype), values=obj.tolist())
+        elif isinstance(obj, np.ndarray) and obj.dtype.kind == 'O':
+            return self.to_dict(obj.tolist())
+        elif isinstance(obj, float):
+            return str(obj) if np.isinf(obj) else obj
+        elif isinstance(obj, int) or isinstance(obj, str) or obj is None:
+            return obj
+        else:
+            raise ValueError("obj %s is not serializable" % str(obj))
+
+
+def to_dict(obj, use_refs=True):
+    return Serializer(use_refs=use_refs).to_dict(obj)
+
+def save(filename, problem):
+    try:
+        p = to_dict(problem)
+        with open(filename, 'w') as fid:
+            json.dump(p, fid)
+    except Exception:
+        traceback.print_exc()
+        warnings.warn(f"failed to create JSON file {filename} for fit problem")
diff --git a/bumps/util.py b/bumps/util.py
index e65b1fc2..2056f152 100644
--- a/bumps/util.py
+++ b/bumps/util.py
@@ -2,12 +2,14 @@
 Miscellaneous utility functions.
 """
 from __future__ import division
+import warnings
 
 __all__ = ["kbhit", "profile", "pushdir", "push_seed", "redirect_console"]
 
 import sys
 import os
 import types
+import inspect
 
 try:  # CRUFT: python 2.x
     from cStringIO import StringIO
@@ -20,6 +22,114 @@
 # so there is no longer a need for bumps.util.erf.
 from scipy.special import erf
 
+# this can be substituted with pydantic dataclass for schema-building...
+try:
+    from typing import Literal, Protocol, runtime_checkable
+except ImportError:
+    from typing_extensions import Literal, Protocol, runtime_checkable
+from typing import Iterable, Optional, Type, TypeVar, Any, Union, Dict, Callable, Tuple, List, Sequence, TYPE_CHECKING
+
+USE_PYDANTIC = os.environ.get('BUMPS_USE_PYDANTIC', "False") == "True"
+if USE_PYDANTIC:
+    from pydantic.dataclasses import dataclass
+    from numpy.typing import NDArray as _NDArray
+    # from pydantic_numpy import NDArray
+else:
+    from dataclasses import dataclass
+    if TYPE_CHECKING:
+        from numpy.typing import NDArray
+
+from dataclasses import field, is_dataclass, Field
+
+def field_desc(description: str) -> Any:
+    return field(metadata={"description": description})
+
+T = TypeVar('T')
+SCHEMA_ATTRIBUTE_NAME = '__bumps_schema__'
+
+def schema(
+        *,
+        include: Optional[List[str]] = None,
+        exclude: Optional[List[str]] = None,
+        classname: Optional[str] = None,
+        eq: bool = True,
+        init: bool = False,
+        frozen: bool = False,
+        repr: bool = True,
+    ) -> Callable[[Type[T]], Type[T]]:
+    
+    """ 
+    Create a dataclass from a subset of field names.
+    fields should be a list of strings corresponding to attribute names,
+    or if fields is None, all annotated attributes will be used as fields
+
+    if attribute "schema_description" is found, will be used for description of
+    generated schema, or else the class docstring will be used.
+    """
+    def set_dataclass(cls: Type[T]) -> Type[T]:
+        realname = cls.__name__
+        if cls.__name__ != cls.__qualname__:
+            warnings.warn(f"serialization does not work on nested classes: {cls.__qualname__} != {cls.__name__}")
+        name = realname if classname is None else classname
+        fqn = f"{cls.__module__}.{name}"
+
+        # check for invalid arguments:
+        # TODO: these arguments appear to never be used... drop them?
+        if include is not None and exclude is not None:
+            raise ValueError(f"{fqn} schema: include array and exclude array are mutually exclusive - only define one")
+
+        dclass_kwargs: Dict[str, Any] = dict(init=init, eq=eq, frozen=frozen, repr=repr)
+
+        # TODO: when pydantic v2 is used, migrate to new function __pydantic_modify_json_schema__
+        # instead of using Config class 
+        # (see https://github.com/pydantic/pydantic/blob/73373c3e08fe5fe23e4b05f549ea34e0da6a16b7/docs/examples/schema_extra_callable.py#L12)
+        if USE_PYDANTIC:
+            # only when generating schema...
+            class Config:
+                @staticmethod
+                def schema_extra(schema: Dict[str, Any], model: Any) -> None:
+                    model.__name__ = name
+                    schema.pop("title")
+                    props = schema.get("properties", {})
+                    if include is not None:
+                        for item in props:
+                            if not item in include:
+                                props.pop(item)
+                    elif exclude is not None:
+                        for item in exclude:
+                            props.pop(item)
+                    props['type'] = {'enum': [fqn]}
+
+            dclass_kwargs["config"] = Config
+            cls.__name__ = name
+
+        setattr(cls, SCHEMA_ATTRIBUTE_NAME, dict(include=include, exclude=exclude, fqn=fqn))
+        dataclass(**dclass_kwargs)(cls)
+        if cls.__name__ != realname:
+            cls.__name__ = realname
+        return cls
+
+    return set_dataclass
+
+def has_schema(cls):
+    return is_dataclass(cls) and hasattr(cls, SCHEMA_ATTRIBUTE_NAME)
+
+@schema(init=True)
+class NumpyArray:
+    """
+    Wrapper for numpy arrays:
+    on serialize, 
+     - array.tolist() is called and stored in 'values' attribute
+     - str(array.dtype) is stored in 'dtype' attribute
+
+     on deserialize,
+      - return new np.ndarray(values, dtype=dtype)
+    """
+    dtype: str
+    values: Sequence = field(default_factory=list)
+
+if USE_PYDANTIC:
+    NDArray = NumpyArray
 
 def parse_errfile(errfile):
     """
@@ -102,6 +212,10 @@ class DynamicModule(types.ModuleType):
     def __init__(self, path, name):
        self.__path__ = [path]
        self.__name__ = name
+       # In the bowels of importlib the parent spec attribute is used to
+       # avoid circular imports, but only if spec is not None. This behaviour
+       # was observed on python 3.11, and perhaps earlier.
+       self.__spec__ = None
 
 
 def relative_import(filename, module_name="relative_import"):
diff --git a/check_fitters.py b/check_fitters.py
index 5578a710..c3fdfc22 100755
--- a/check_fitters.py
+++ b/check_fitters.py
@@ -74,11 +74,11 @@ def check_fit(fitter, store, targets):
     model_index = 0
     with open(errfiles[0]) as fid:
         for line in fid:
-            if line.startswith("[chisq="):
-                if line[7:10].lower() == 'inf':
+            if line.startswith("[overall chisq="):
+                if line[15:10].lower() == 'inf':
                     value = np.inf
                 else:
-                    value = float(line[7:].split("(")[0])
+                    value = float(line[15:].split("(")[0])
                 assert abs(value-targets[model_index]) < 1e-2, \
                     "error in %s: expected %.3f but got %.3f" \
                     % (fitter, targets[model_index], value)
diff --git a/doc/conf.py b/doc/conf.py
index aae12795..2a18d44b 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -62,6 +62,7 @@
 
 nitpick_ignore = [
     ('py:class', 'type'),
+    ('py:class', 'enum.Enum'),
     ('py:class', 'object'),
     ("py:class", "numpy.ndarray"),
     ("py:class", "np.ndarray"),
diff --git a/doc/guide/intro.rst b/doc/guide/intro.rst
index 7db1df23..11f8a15b 100644
--- a/doc/guide/intro.rst
+++ b/doc/guide/intro.rst
@@ -39,14 +39,14 @@ the line.
 Every fit file ends with a :func:`FitProblem <bumps.fitproblem.FitProblem>`
 definition::
 
-    problem = FitProblem(model)
+    problem = FitProblem(models)
 
 In fact, this is the only requirement of the fit file.  The Bumps engine
 loads the fit file, retieves the *problem* symbol and feeds it to one of the
 :mod:`fitters <bumps.fitters>`.  Some fit files do not even use *FitProblem* to
 define *problem*, or use *Parameter* objects for the fitted parameters, so
 long as *problem* implements the
-:class:`BaseFitProblem <bumps.fitproblem.BaseFitProblem>` interface, which
+:class:`FitProblem <bumps.fitproblem.FitProblem>` interface, which
 provides *getp* to get the existing parameter vector, *setp* to set a new
 parameter vector, *bounds* to return the parameter bounds, and *nllf* to
 to compute the negative log likelihood function.  The remaining methods are
@@ -60,7 +60,7 @@ pattern to use is::
 
     import bumps.names as bmp
     ...
-    problem = bmp.FitProblem(model)
+    problem = bmp.FitProblem(models)
 
 This documents to the reader unfamiliar with your code (such as you, dear
 reader, when looking at your model files two years from now) exactly where
diff --git a/doc/guide/parameter.rst b/doc/guide/parameter.rst
index 2142a055..262b2411 100644
--- a/doc/guide/parameter.rst
+++ b/doc/guide/parameter.rst
@@ -76,8 +76,6 @@ are helpful:
 - :meth:`soft_range <Parameter.soft_range>` is a combination of *range*
   and *dev* in that the parameter has equal probability within [*low*,*high*]
   but Gaussian probability of width *std* as it strays outside of the range.
-- :meth:`pdf <Parameter.pdf>` is like *dev* but works with any continuous
-  `scipy.stats distribution <http://docs.scipy.org/doc/scipy/reference/stats.html>`_.
 
 All these methods set the *bounds* attribute on the parameter in one way
 or another.  See :mod:`bumps.bounds` for details.  Technically, setting
diff --git a/extra/sasview/sasbumps.py b/extra/sasview/sasbumps.py
index 79d44f55..16f30d71 100644
--- a/extra/sasview/sasbumps.py
+++ b/extra/sasview/sasbumps.py
@@ -7,7 +7,7 @@
 import sys
 import numpy as np
 from bumps.names import Parameter, FitProblem, FreeVariables, pmath, preview, fit
-from bumps.parameter import BaseParameter
+from bumps.parameter import Parameter
 from bumps.util import push_seed as seed
 import sans
 import sans.models
@@ -52,7 +52,7 @@ def load_model(model, name=None, **kw):
                            %(k, model, formatted_pars))
         if '.' in k and not k.endswith('.width'):
             M.setParam(k, v)
-        elif isinstance(v, BaseParameter):
+        elif isinstance(v, Parameter):
             M._bumps_pars[k] = v
         elif isinstance(v, (tuple, list)):
             low, high = v
diff --git a/setup.py b/setup.py
index 4a81ea3a..183b47f0 100755
--- a/setup.py
+++ b/setup.py
@@ -51,7 +51,8 @@
     packages=packages,
     package_data=gui_resources.package_data(),
     scripts=scripts,
-    install_requires=['six'],
+    # TODO: typing_extensions is only required for python 3.7 and below
+    install_requires=['typing_extensions>=3.7.4', 'dataclasses; python_version < "3.7.0"'],
     #install_requires=['six', 'numpy>=1.0', 'scipy>=0.7.0', 'matplotlib>=1.0'],
     #install_requires = ['httplib2'],
 )


Test patch:
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 0e98201f..49641985 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -10,8 +10,9 @@ jobs:
       matrix:
         cfg:
           #- { os: ubuntu-latest, py: 2.7 }
-          - { os: ubuntu-latest, py: 3.6 }
           - { os: ubuntu-latest, py: 3.8, doc: 1 }
+          - { os: ubuntu-latest, py: "3.10" }
+          - { os: ubuntu-latest, py: 3.11 }
           - { os: windows-latest, py: 3.8 }
           - { os: macos-latest, py: 3.8 }
       fail-fast: false
@@ -34,6 +35,11 @@ jobs:
         python -m pip install --upgrade pip
         python -m pip install wheel setuptools
         python -m pip install numpy scipy matplotlib scikit-learn pytest pytest-cov
+    
+    - name: Install dataclasses for python3.6
+      if: matrix.cfg.py == 3.6
+      run: |
+        python -m pip install dataclasses
 
     - name: Build package [Linux only]
       if: ${{ runner.os == 'Linux' }}
diff --git a/pytest.ini b/pytest.ini
index be4e4cc7..94bc40a6 100644
--- a/pytest.ini
+++ b/pytest.ini
@@ -2,7 +2,7 @@
 addopts = --doctest-modules --doctest-glob=*.rst --cov=bumps
 doctest_optionflags = ELLIPSIS NORMALIZE_WHITESPACE
 testpaths = bumps doc/getting_started doc/guide
-norecursedirs = numdifftools gui
+norecursedirs = numdifftools gui webview/client
 python_files = *.py
 python_classes = NoClassTestsWillMatch
-python_functions = test *_test test_*
\ No newline at end of file
+python_functions = test *_test test_*