numerical

DummyNumExpr

Bases: NumExprBase

Dummy callable object which returns supplied 'obj' when called Has no parameters or symbols

Source code in slimfit/numerical.py

class DummyNumExpr(NumExprBase):
    """Dummy callable object which returns supplied 'obj' when called
    Has no parameters or symbols
    """

    def __init__(self, obj: Any, *args, **kwargs):
        #
        super().__init__(*args, **kwargs)
        self.obj = obj

    def __call__(self, **kwargs):
        return self.obj

    @property
    def symbols(self) -> set[Symbol]:
        return set()

    @property
    def shape(self) -> Shape:
        try:
            return self.obj.shape
        except AttributeError:
            return ()

GMM

Bases: CompositeExpr

Source code in slimfit/numerical.py

class GMM(CompositeExpr):
    # todo can also be implemented as normal NumExpr but with broadcasting parameter shapes
    # important is that GMM class allows users to find oud positions of parmaeters in mu / sigma
    # matrices for EM GMM optimization.

    def __init__(
        self,
        x: Symbol | NumExpr,
        mu: Matrix | MatrixNumExpr,
        sigma: Matrix | MatrixNumExpr,
        name: Optional[str] = None,
    ):
        # check for the correct shape when creating from sympy Matrix
        if isinstance(mu, Matrix) and mu.shape[0] != 1:
            raise ValueError(
                "GMM parameter matrices must be of shape (1, N) where N is the number of states."
            )

        expr = {"x": x, "mu": mu, "sigma": sigma}

        # todo some kind of properties object for this metadata
        name = name or "GMM"  # counter for number of instances?
        self.kind = "gmm"
        super().__init__(expr)

    def __call__(self, **kwargs) -> np.ndarray:
        result = super().__call__(**kwargs)

        x, mu, sig = result["x"], result["mu"], result["sigma"]
        output = 1 / (np.sqrt(2 * np.pi) * sig) * np.exp(-np.power((x - mu) / sig, 2) / 2)

        # Output shape is (datapoints, states, 1) to match output shape of matrix
        # exponentiation model
        return np.expand_dims(output, -1)

    def __repr__(self):
        return f"GMM({self.expr['x']}, {self.expr['mu']}, {self.expr['sigma']})"

    def to_numerical(self) -> GMM:
        # todo probably this is the same as the super class method
        num_expr = {k: to_numerical(expr) for k, expr in self.items()}
        instance = GMM(**num_expr)

        return instance

    @property
    def shape(self) -> Shape:
        shape = super().shape
        return shape + (1,)

    @property
    def states(self) -> Optional[list[str]]:
        """
        List of state names, if naming scheme of mu's is of format `mu_<state_name>`.
        """

        mus, suffices = zip(*(elem.name.split("_") for elem in self["mu"]))

        if all((mu == "mu" for mu in mus)):
            return list(suffices)
        else:
            return None

states property

List of state names, if naming scheme of mu's is of format mu_<state_name>.

MarkovIVP

Bases: CompositeExpr

Uses scipy.integrate.solve_ivp to numerically find time evolution of a markov process given a transition rate matrix.

Returned shape is (len(t_var), len(y0), 1), or (, , 1)

Source code in slimfit/numerical.py

class MarkovIVP(CompositeExpr):
    """Uses scipy.integrate.solve_ivp to numerically find time evolution of a markov process
        given a transition rate matrix.

    Returned shape is (len(t_var), len(y0), 1), or (<datapoints>, <states>, 1)

    """

    def __init__(
        self,
        t: Symbol | NumExpr,
        trs_matrix: Matrix | MatrixNumExpr,
        y0: Matrix | MatrixNumExpr,
        domain: Optional[tuple[float, float]] = None,
        **ivp_kwargs,
    ):
        expr = {"t": t, "trs_matrix": trs_matrix, "y0": y0}

        super().__init__(expr)

        ivp_defaults = {"method": "Radau"}
        self.ivp_defaults = ivp_defaults | ivp_kwargs
        self.domain = domain

    def __call__(self, **kwargs):
        result = super().__call__(**kwargs)

        # if `self['t']` does not depend on any parameters; domain can be precomputed and
        # does not have to be determined for every call
        # although its every fast to do so
        domain = self.domain or self.get_domain(result["t"])
        sol = solve_ivp(
            self.grad_func,
            domain,
            y0=result["y0"].squeeze(),
            t_eval=result["t"],
            args=(result["trs_matrix"],),
            **self.ivp_defaults,
        )

        # shape is modified to match the output shape of matrix exponentiation model
        # exp(m*t) @ y0, which is (datapoints, states, 1)
        return np.expand_dims(sol.y.T, -1)

    # TODO generalize
    def to_numerical(self) -> MarkovIVP:
        num_expr = {k: to_numerical(expr) for k, expr in self.items()}
        instance = MarkovIVP(**num_expr, domain=self.domain, **self.ivp_defaults)

        return instance

    def get_domain(self, arr: np.ndarray) -> tuple[float, float]:
        # padding?
        return arr[0], arr[-1]

    @staticmethod
    def grad_func(t, y, trs_matrix):
        return trs_matrix @ y

MatrixNumExpr

Bases: NumExpr

Source code in slimfit/numerical.py

class MatrixNumExpr(NumExpr):
    def __init__(
        self,
        expr: MatrixBase,
        name: Optional[str] = None,
        kind: Optional[str] = None,
    ):
        if not isinstance(expr, MatrixBase):
            raise TypeError("Expression must be an instance of MatrixParameter or sympy.Matrix")
        self._name = name

        # Callable type is used by minimizers to determine solving strategy
        if kind is None:
            self.kind = identify_expression_kind(expr)
        elif isinstance(kind, str):
            self.kind: str = kind.lower()
        else:
            raise TypeError("Invalid type for 'kind', must be 'str'")

        super().__init__(expr)

    @property
    def name(self) -> str:
        if self._name:
            return self._name
        if self.kind == "constant":
            symbol_names = [str(symbol) for symbol in self.expr]
            prefix = [name.split("_")[0] for name in symbol_names]
            if len(set(prefix)) == 1:  # All prefixes are identical, prefix is the name
                return prefix[0]
        else:
            return "M"

    @name.setter
    def name(self, name):
        self._name = name

    @property
    def shape(self) -> Shape:
        # again there might be problems depending on how the matrix elements depend on
        # different combinations of parameters and data
        # for now we assume this is the same for all elements

        # Find the shape from broadcasting parameters and data
        base_shape = super().shape

        # squeeze last dim if shape is (1,)
        base_shape = () if base_shape == (1,) else base_shape
        shape = base_shape + self.expr.shape

        return shape

    @cached_property
    def element_mapping(self) -> dict[Symbol, tuple[int, int]]:
        """
        Dictionary mapping Symbol to matrix indices
        only returns entries if the matrix element is equal to a `Symbol`, not `Expr`
        """

        element_mapping = {}
        for i, j in np.ndindex(self.expr.shape):
            elem = self.expr[i, j]
            if isinstance(elem, Symbol):
                element_mapping[elem] = (i, j)
        return element_mapping

    @cached_property
    def lambdified(self) -> np.ndarray:
        """Array of lambdified function per matrix element"""
        # TODO scalercallable per element

        lambdas = np.empty(self.expr.shape, dtype=object)
        # todo might go wrong when not all elements have the same parameters
        for i, j in np.ndindex(self.expr.shape):
            lambdas[i, j] = lambdify(sorted(self.symbols, key=str), self.expr[i, j])

        return lambdas

    def forward(self, **kwargs):
        """forward pass with type checking

        or without?
        """
        ...

    def __call__(self, **kwargs: float) -> np.ndarray:
        # https://github.com/sympy/sympy/issues/5642
        # Prepare kwargs for lambdified
        # try:
        #     parameters: dict[str, np.ndarray | float] = {
        #         k: kwargs[k] for k in self.free_parameters.keys()
        #     }
        # except KeyError as e:
        #     raise KeyError(f"Missing value for parameter {e}") from e

        # check shapes
        # this should move somewhere else
        # for p_name, p_value in parameters.items():
        #     if getattr(p_value, "shape", tuple()) != self.parameters[p_name].shape:
        #         raise ValueError(f"Shape mismatch for parameter {p_name}")

        ld_kwargs = self.parse_kwargs(**kwargs)

        # todo precomputed shapes
        # try:
        #     # Shape is given be pre-specified shape
        #     shape = self.shape
        # except AttributeError:
        base_shape = np.broadcast_shapes(
            *(getattr(value, "shape", tuple()) for value in ld_kwargs.values())
        )

        # squeeze last dim if shape is (1,)
        base_shape = () if base_shape == (1,) else base_shape
        shape = base_shape + self.expr.shape

        out = np.empty(shape)
        for i, j in np.ndindex(self.expr.shape):
            out[..., i, j] = self.lambdified[i, j](**ld_kwargs)

        return out

    @property
    def values(self) -> np.ndarray:
        """

        Returns: Array with elements set to parameter values

        """
        raise DeprecationWarning("Deprecate in favour of `value`")
        arr = np.empty(self.shape)
        for i, j in np.ndindex(self.shape):
            arr[i, j] = self.expr[i, j].value

        return arr

    def __getitem__(self, key):
        return self.expr[key]

    def __contains__(self, item) -> bool:
        return self.expr.__contains__(item)

    # this is more of a str than a repr
    def __repr__(self):
        names = sorted(self.symbol_names)
        name = self.name or "MatrixNumExpr"
        return f"{name}({', '.join(names)})"

    def index(self, symbol: Symbol) -> tuple[int, int]:
        """
        Returns indices of parameter for Matrix Expressions

        Args:
            name: Parameter name to find matrix elements of

        Returns: Tuple of matrix elements ij

        """

        return self.element_mapping[symbol]

element_mapping cached property

Dictionary mapping Symbol to matrix indices only returns entries if the matrix element is equal to a Symbol, not Expr

lambdified cached property

Array of lambdified function per matrix element

values property

Returns: Array with elements set to parameter values

forward(**kwargs)

forward pass with type checking

or without?

Source code in slimfit/numerical.py

def forward(self, **kwargs):
    """forward pass with type checking

    or without?
    """
    ...

index(symbol)

Returns indices of parameter for Matrix Expressions

Parameters:

Name	Type	Description	Default
name		Parameter name to find matrix elements of	required

Returns: Tuple of matrix elements ij

Source code in slimfit/numerical.py

def index(self, symbol: Symbol) -> tuple[int, int]:
    """
    Returns indices of parameter for Matrix Expressions

    Args:
        name: Parameter name to find matrix elements of

    Returns: Tuple of matrix elements ij

    """

    return self.element_mapping[symbol]

NumberExpr

Bases: NumExprBase

callable object for sympy Numbers

Source code in slimfit/numerical.py

class NumberExpr(NumExprBase):
    """callable object for sympy `Number`s"""

    def __init__(self, obj: Any, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.obj = obj

    def __call__(self, **kwargs):
        return self.obj.evalf()

    @property
    def symbols(self) -> set[Symbol]:
        return set()

    @property
    def shape(self) -> Shape:
        try:
            return self.obj.shape
        except AttributeError:
            return ()

identify_expression_kind(sympy_expression)

Find the type of expression

Only implemented for 'constant' kind, otherwise return is generic

Source code in slimfit/numerical.py

def identify_expression_kind(sympy_expression: Union[Expr, MatrixBase]) -> str:
    """Find the type of expression

    Only implemented for 'constant' kind, otherwise return is generic

    """

    # check for gaussian mixture model ...

    if isinstance(sympy_expression, MatrixBase):
        ...

        if all(isinstance(elem, Symbol) for elem in sympy_expression):
            # this should also check the symbols for their shapes, for it to be constant
            # the elements should all be scalars
            return "constant"

    return "generic"

to_numerical(expression)

Converts sympy expression to slimfit numerical expression

if the expressions already is an NumExpr; the object is modified in-place by setting the parameters

Source code in slimfit/numerical.py

def to_numerical(
    expression: Union[NumExprBase, Expr, MatrixBase | Model | CompositeExpr],
) -> NumExprBase | CompositeExpr:
    """Converts sympy expression to slimfit numerical expression

    if the expressions already is an NumExpr; the object is modified in-place by setting
    the parameters

    """

    if hasattr(expression, "to_numerical"):
        return expression.to_numerical()
    # elif isinstance(expression, Model):
    #     model_dict = {lhs: to_numerical(rhs) for lhs, rhs in expression.items()}
    #     from slimfit.models import NumericalModel
    #
    #     return NumericalModel(model_dict, parameters, data)
    if isinstance(expression, HadamardProduct):
        from slimfit.operations import Mul

        return Mul(*(to_numerical(arg) for arg in expression.args))
    elif isinstance(expression, MatrixBase):
        return MatrixNumExpr(expression)
    elif isinstance(expression, Expr):
        return NumExpr(expression)
    elif isinstance(expression, (np.ndarray, float, int)):
        return DummyNumExpr(expression)
    elif isinstance(expression, NumExprBase):
        return expression
    else:
        raise TypeError(f"Invalid type {type(expression)!r}")