# Bumps 1.0.0a13 documentation

## bumps.parameter

# Source code for bumps.parameter

    # This program is public domain
    # Author: Paul Kienzle
    """
    Fitting parameter objects.

    Parameters are a big part of the interface between the model and the fitting
    engine.  By saving and retrieving values and ranges from the parameter, the
    fitting engine does not need to be aware of the structure of the model.

    Users can also perform calculations with parameters, tying together different
    parts of the model, or different models.
    """

    # __all__ = [ 'Parameter']
    import operator
    import sys
    import builtins
    from dataclasses import dataclass, field, InitVar
    from functools import reduce
    import warnings
    from copy import copy
    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_desc, schema_config

    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
    # calculation via a parameterized function.  This will require a restructuring
    # of the parameter claas.  The park-1.3 solution is viable: given a parameter
    # set, figure out which order the expressions need to be evaluated by
    # building up a dependency graph.  With a little care, we can check which
    # parameters have actually changed since the last calculation update, and
    # restrict the dependency graph to just them.
    # TODO: support full aliasing, so that floating point model attributes can
    # be aliased to a parameter.  The same technique as subexpressions applies:
    # when the parameter is changed, the model will be updated and will need
    # to be re-evaluated.

    # TODO: maybe move this to util?

    def to_dict(p):
        if hasattr(p, "to_dict"):
            return p.to_dict()
        elif isinstance(p, (tuple, list)):
            return [to_dict(v) for v in p]
        elif isinstance(p, dict):
            return {k: to_dict(v) for k, v in p.items()}
        elif isinstance(p, (bool, str, float, int, type(None))):
            return p
        elif isinstance(p, np.ndarray):
            # TODO: what about inf, nan and object arrays?
            return p.tolist()
        elif False and callable(p):
            # TODO: consider including functions and arbitrary values
            import base64
            import dill

            encoding = base64.encodebytes(dill.dumps(p)).decode("ascii")
            return {"type": "dill", "value": str(p), "encoding": encoding}
            ## To recovert the function
            # if allow_unsafe_code:
            #     encoding = item['encoding']
            #     p = dill.loads(base64.decodebytes(encoding).encode('ascii'))
        else:
            # print(f"converting type {type(p)} to str")
            return str(p)

    @dataclass(init=False)
    class Uniform:
        """Uniform distribution with hard boundaries"""

    @dataclass(init=False)
    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
    # @dataclass(init=False)

    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:
        """
        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

    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 []

    @dataclass(init=False)
    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_config()
    @dataclass(init=False)
    class Parameter(ValueProtocol, 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]
        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
        tags: List[str] = field(default_factory=list)

        _fixed: bool

        def parameters(self):
            pars = [self]
            if hasattr(self.slot, "parameters"):
                pars += self.slot.parameters()
            return pars

        def pmp(self, plus, minus=None, limits=None):
            """
            Allow the parameter to vary as value +/- percent.

            pmp(*percent*) -> [value*(1-percent/100), value*(1+percent/100)]

            pmp(*plus*, *minus*) -> [value*(1+minus/100), value*(1+plus/100)]

            In the *plus/minus* form, one of the numbers should be plus and the
            other minus, but it doesn't matter which.

            If *limits* are provided, bound the end points of the range to lie
            within the limits.

            The resulting range is converted to "nice" numbers.
            """
            bounds = mbounds.pmp(self.value, plus, minus, limits=limits)
            self.bounds = bounds
            self.fixed = False
            return self

        def pm(self, plus, minus=None, limits=None):
            """
            Allow the parameter to vary as value +/- delta.

            pm(*delta*) -> [value-delta, value+delta]

            pm(*plus*, *minus*) -> [value+minus, value+plus]

            In the *plus/minus* form, one of the numbers should be plus and the
            other minus, but it doesn't matter which.

            If *limits* are provided, bound the end points of the range to lie
            within the limits.

            The resulting range is converted to "nice" numbers.
            """
            bounds = mbounds.pm(self.value, plus, minus, limits=limits)
            self.bounds = bounds
            self.fixed = False
            return self

        def dev(self, std, mean=None, limits=None, sigma=None, mu=None):
            """
            Allow the parameter to vary according to a normal distribution, with
            deviations from the mean added to the overall cost function for the
            model.

            If *mean* is None, then it defaults to the current parameter value.

            If *limits* are provide, then use a truncated normal distribution.

            Note: *sigma* and *mu* have been replaced by *std* and *mean*, but
            are left in for backward compatibility.
            """
            if sigma is not None or mu is not None:
                # CRUFT: remove sigma and mu parameters
                warnings.warn(DeprecationWarning("use std,mean instead of mu,sigma in Parameter.dev"))
                if sigma is not None:
                    std = sigma
                if mu is not None:
                    mean = mu
            if mean is None:
                mean = self.value  # Note: value is an attribute of the derived class
            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.
        #     """
        #     # 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 = (low, high)
            self.distribution = Uniform()
            self.fixed = False
            return self

        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 = (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 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 = 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

        def __float__(self):
            return float(self.value)

        def nllf(self) -> float:
            """
            Return -log(P) for the current parameter 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) -> float:
            """
            Return the z score equivalent for the current parameter value.

            That is, the given the value of the parameter in the underlying
            distribution, find the equivalent value in the standard normal.
            For a gaussian, this is the z score, in which you subtract the
            mean and divide by the standard deviation to get the number of
            sigmas away from the mean.  For other distributions, you need to
            compute the cdf of value in the parameter distribution and invert
            it using the ppf from the standard normal distribution.
            """
            return 0.0 if self.prior is None else self.prior.residual(self.value)

        def valid(self):
            """
            Return true if the parameter is within the valid range.
            """
            return not isinf(self.nllf())

        def format(self):
            """
            Format the parameter, value and range as a string.
            """
            return "%s=%g in %s" % (self, self.value, self.prior)

        def __str__(self):
            name = self.name if self.name is not None else "?"
            return name

        def __repr__(self):
            return "Parameter(%s)" % self

        # TODO: deprecate

        @classmethod
        def default(cls: type, value: Union[float, Tuple[float, float], ValueType], **kw) -> "Parameter":
            """
            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, ValueProtocol):
                return value
            else:
                return cls(value, **kw)

        def set(self, value):
            """
            Set a new value for the parameter, ignoring the bounds.
            """
            self.slot.value = value

        def clip_set(self, value):
            """
            Set a new value for the parameter, clipping it to the bounds.
            """
            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,
            tags: Optional[List[str]] = None,
            **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
                    value = None
                except TypeError:
                    pass
            if fixed is None:
                fixed = bounds is None
            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))

            self.slot = slot
            self.name = name
            self.id = id if id is not None else str(uuid.uuid4())
            self.tags = tags if tags is not None else []
            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 = (
                    (-np.inf if bounds[0] is None else float(bounds[0])),
                    (np.inf if bounds[1] is None else float(bounds[1])),
                )
            self.bounds = bounds
            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.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.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.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 add_tag(self, tag: str):
            if not tag in self.tags:
                self.tags.append(tag)

        def remove_tag(self, tag: Optional[str] = None):
            if tag is None:
                self.tags = []
            else:
                self.tags = [t for t in self.tags if not t == tag]

        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

    def tag_all(parameter_tree, tag, remove=False):
        if isinstance(parameter_tree, dict):
            tag_all([item for item in parameter_tree.values()], tag, remove=remove)
        elif hasattr(parameter_tree, "add_tag"):
            if remove:
                parameter_tree.remove_tag(tag)
            else:
                parameter_tree.add_tag(tag)
        elif hasattr(parameter_tree, "parameters"):
            tag_all(parameter_tree.parameters(), tag, remove=remove)
        elif hasattr(parameter_tree, "__iter__"):
            for item in parameter_tree:
                tag_all(item, tag, remove=remove)
        else:
            warnings.warn(f"parameter tree should have only list, object and Parameter items: {parameter_tree}")

    def untag_all(parameter_tree, tag: Optional[str] = None):
        tag_all(parameter_tree, tag, remove=True)

    @dataclass
    class Variable(ValueProtocol):
        """
        Saved state for a random variable in the model.
        """

        value: float

        def parameters(self):
            return []

    @schema_config()
    @dataclass(init=True, frozen=True, eq=False)
    class Constant(ValueProtocol):  # type: ignore
        """
        Saved state for an unmodifiable value.

        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.
        """

        value: float
        name: Optional[str] = None
        id: str = field(metadata={"format": "uuid"}, default_factory=lambda: str(uuid.uuid4()))

        fittable = False  # class property fixed across all objects
        fixed = True  # class property fixed across all objects

        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 = "truediv"
        floordiv = "floordiv"
        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,
        "truediv": 6,
        "floordiv": 6,
        "add": 7,
        "sub": 7,
        "gt": 12,
        "lt": 12,
        "ge": 12,
        "le": 12,
        "eq": 12,
        "ne": 12,
    }
    OPERATOR_STRING = {
        "pow": "**",
        "pos": "+",
        "neg": "-",
        "mul": "*",
        "truediv": "/",
        "floordiv": "//",
        "add": "+",
        "sub": "-",
        "gt": ">",
        "lt": "<",
        "ge": ">=",
        "le": "<=",
        "eq": "==",
        "ne": "!=",
    }

    def _lookup_operator(op_name):
        if not hasattr(Operators, op_name) and op_name not in UserFunctionRegistry:
            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 UserFunctionRegistry
            fn = UserFunctionRegistry.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.

        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

    @dataclass(init=False)
    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 self._fn(*(float(arg) for arg in self.args))

        @property
        def name(self):
            return str(self)

        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", "floordiv", "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()

    UserFunctionRegistry: Dict[str, Callable[..., float]] = {}

    # TODO: allow schema validation on user-defined functions

    @dataclass(init=False)
    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

        # 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.
        def __init__(self, fn: Callable):
            name = fn.__name__
            if name in UserFunctionRegistry:
                raise TypeError(f"Function {name} already registered in bumps.")
            UserFunctionRegistry[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))

    @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]

    @dataclass(init=False)
    class ParameterSet:
        """
        A parameter that depends on the model.
        """

        names: Optional[List[str]]
        reference: Parameter
        parameterlist: Optional[List[Parameter]]

        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.

            *names* is the list of model names.

            *reference* is the underlying :class:`parameter.Parameter` that will
            be set when the model is selected.

            *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
            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))

            # 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 {
                "type": "ParameterSet",
                "names": self.names,
                "reference": to_dict(self.reference),
                # Note: parameters are stored in a numpy array
                "parameters": to_dict(self.parameters.tolist()),
            }

        # Make the parameter set act like a list
        def __getitem__(self, i):
            """
            Return the underlying parameter for the model index.  Index can
            either be an integer or a model name.  It can also be a slice,
            in which case a new parameter set is returned.
            """
            # Try looking up the free variable by model name rather than model
            # index. If this fails, assume index is a model index.
            try:
                i = self.names.index(i)
            except ValueError:
                pass
            if isinstance(i, slice):
                obj = copy(self)
                obj.names = self.names[i]
                obj.reference = self.reference
                obj.parameters = self.parameters[i]
                return obj
            return self.parameters[i]

        def __setitem__(self, i, v):
            """
            Set the underlying parameter for the model index.  Index can
            either be an integer or a model name.  It can also be a slice,
            in which case all underlying parameters are set, either to the
            same value if *v* is a single parameter, otherwise *v* must have
            the same length as the slice.
            """
            try:
                i = self.names.index(i)
            except ValueError:
                pass
            self.parameters[i] = v

        def __iter__(self):
            return iter(self.parameters)

        def __len__(self):
            return len(self.parameters)

        def set_model(self, index):
            """
            Set the underlying model parameter to the value of the nth model.
            """
            self.reference.value = self.parameters[index].value

        def get_model(self, index):
            """
            Get the reference and underlying model parameter for the nth model.
            """
            return (id(self.reference), self.parameters[index])

        @property
        def values(self):
            return [p.value for p in self.parameters]

        @values.setter
        def values(self, values):
            for p, v in zip(self.parameters, values):
                p.value = v

        def range(self, *args, **kw):
            """
            Like :meth:`Parameter.range`, but applied to all models.
            """
            for p in self.parameters:
                p.range(*args, **kw)

        def pm(self, *args, **kw):
            """
            Like :meth:`Parameter.pm`, but applied to all models.
            """
            for p in self.parameters:
                p.pm(*args, **kw)

        def pmp(self, *args, **kw):
            """
            Like :meth:`Parameter.pmp`, but applied to all models.
            """
            for p in self.parameters:
                p.pmp(*args, **kw)

    class Reference(Parameter):
        """
        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.

        NOTE: this class can not be serialized with a dataclass schema
        TODO: can sasmodels just use Parameter directly?
        """

        def __init__(self, obj, attr, **kw):
            self.obj = obj
            self.attr = attr
            kw.setdefault("name", ".".join([obj.__class__.__name__, attr]))
            Parameter.__init__(self, **kw)

        @property
        def value(self):
            return getattr(self.obj, self.attr)

        @value.setter
        def value(self, value):
            setattr(self.obj, self.attr, value)

    @dataclass(init=False)
    class FreeVariables(object):
        """
        A collection of parameter sets for a group of models.

        *names* is the set of model names.

        The parameters themselves are specified as key=value pairs, with key
        being the attribute name which is used to retrieve the parameter set
        and value being a :class:`Parameter` containing the parameter that is
        shared between the models.

        In order to evaluate the log likelihood of all models simultaneously,
        the fitting program will need to call set_model with the model index
        for each model in turn in order to substitute the values from the free
        variables into the model.  This allows us to share a common sample
        across multiple data sets, with each dataset having its own values for
        some of the sample parameters.  The alternative is to copy the entire
        sample structure, sharing references to common parameters and creating
        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, parametersets=None, **kw):
            if names is None:
                raise TypeError("FreeVariables needs name=[model1, model2, ...]")
            self.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
        # chokes without it.
        def __getstate__(self):
            return self.__dict__

        def __setstate__(self, state):
            self.__dict__ = state

        def __getattr__(self, k):
            """
            Return the parameter set for the given free parameter.
            """
            try:
                return self.parametersets[k]
            except KeyError:
                raise AttributeError("FreeVariables has no attribute %r" % k)

        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())

        def to_dict(self):
            return {"type": type(self).__name__, "names": self.names, "parameters": to_dict(self.parametersets)}

        def set_model(self, i):
            """
            Set the reference parameters for model *i*.
            """
            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())

    def flatten(s):
        if isinstance(s, (tuple, list, np.ndarray)):
            return reduce(lambda a, b: a + flatten(b), s, [])
        elif isinstance(s, set):
            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, ValueProtocol):
            return [s]
        elif s is None:
            return []
        else:
            raise TypeError("don't understand type %s for %r" % (type(s), s))

    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 = []
            for k in sorted(p.keys()):
                if k.startswith("_"):
                    continue
                s = format(p[k], indent + 2, field=k, freevars=freevars)
                label = " " * indent + "." + k
                if s.endswith("\n"):
                    res.append(label + "\n" + s)
                else:
                    res.append(label + " = " + s + "\n")
            if "_index" in p:
                res.append(format(p["_index"], indent, freevars=freevars))
            return "".join(res)

        elif isinstance(p, (list, tuple, np.ndarray)) and len(p):
            res = []
            for k, v in enumerate(p):
                s = format(v, indent + 2, freevars=freevars)
                label = " " * indent + "[%d]" % k
                if s.endswith("\n"):
                    res.append(label + "\n" + s)
                else:
                    res.append(label + " = " + s + "\n")
            return "".join(res)

        elif isinstance(p, Parameter):
            s = ""
            if str(p) != field:
                s += str(p) + " = "
            s += "%g" % p.value
            if not p.fixed:
                if p.prior is not None:
                    bounds = p.prior.limits
                elif p.bounds is not None:
                    bounds = p.bounds
                else:
                    bounds = p.limits
                s += " in [%g,%g]" % tuple(bounds)
            return s

        elif isinstance(p, Parameter):
            return "%s = %g" % (str(p), p.value)

        else:
            return str(p)

    def summarize(pars, sorted=False):
        """
        Return a stylized list of parameter names and values with range bars
        suitable for printing.

        If sorted, then print the parameters sorted alphabetically by name.
        """
        output = []
        if sorted:
            pars = sorted(pars, key=lambda x: x.name)
        for p in pars:
            if not isfinite(p.value):
                bar = ["*invalid* "]
            else:
                position = int(p.prior.get01(p.value) * 9.999999999)
                bar = ["."] * 10
                if position < 0:
                    bar[0] = "<"
                elif position > 9:
                    bar[9] = ">"
                else:
                    bar[position] = "|"
            output.append("%40s %s %10g in %s" % (p.name, "".join(bar), p.value, p.bounds))
        return "\n".join(output)

    def unique(s) -> List[Parameter]:
        """
        Return the unique set of parameters

        The ordering is stable.  The same parameters/dependencies will always
        return the same ordering, with the first occurrence first.
        """
        # Walk structures such as dicts and lists
        pars = flatten(s)
        # print "====== flattened"
        # print "\n".join("%s:%s"%(id(p),p) for p in pars)
        # Also walk parameter expressions
        pars = pars + flatten([p.parameters() for p in pars])
        # print "====== extended"
        # print "\n".join("%s:%s"%(id(p),p) for p in pars)

        # TODO: implement n log n rather than n^2 uniqueness algorithm
        # problem is that the sorting has to be unique across a pickle.
        result = []
        for p in pars:
            if not any(p is q for q in result):
                result.append(p)

        # print "====== unique"
        # print "\n".join("%s:%s"%(id(p),p) for p in result)
        # Return the complete set of parameters
        return result

    def fittable(s):
        """
        Return the list of fittable parameters in no paraticular order.

        Note that some fittable parameters may be fixed during the fit.
        """
        return [p for p in unique(s) if p.fittable]

    def varying(s: List[Parameter]) -> List[Parameter]:
        """
        Return the list of fitted parameters in the model.

        This is the set of parameters that will vary during the fit.
        """
        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: 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.prior.random(1)[0]

    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 = '!='

    @dataclass(init=False)
    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.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 _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):
        """
        Parameter alias.

        Rather than modifying a model to contain a parameter slot,
        allow the parameter to exist outside the model. The resulting
        parameter will have the full parameter semantics, including
        the ability to replace a fixed value with a parameter expression.

        """

        def __init__(self, obj, attr, p=None, name=None):
            self.obj = obj
            self.attr = attr
            if name is None:
                name = ".".join([obj.__class__.__name__, attr])
            self.p = Parameter.default(p, name=name)

        def update(self):
            setattr(self.obj, self.attr, self.p.value)

        def parameters(self):
            return self.p.parameters()

        def to_dict(self):
            return {
                "type": type(self).__name__,
                "p": to_dict(self.p),
                # TODO: can't json pickle arbitrary objects
                "obj": to_dict(self.obj),
                "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")
        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.0
        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()