minimizers

ConstantOptimizer

Bases: EMOptimizer

model is of form {Probability(px): Matrix[[]] ...} where all matrix elements are just scalar parameters.

Source code in slimfit/minimizers.py

class ConstantOptimizer(EMOptimizer):

    """
    model is of form {Probability(px): Matrix[[<parameters>]] ...} where all matrix elements are just scalar parameters.
    """

    def step(self) -> dict[str, float]:
        parameters = {}
        # for p_name in self.model.free_parameters:
        for parameter in self.free_parameters:
            num, denom = 0.0, 0.0
            for lhs, rhs in self.model.items():
                # rhs is of type MatrixNumExpr
                if parameter.symbol in rhs.symbols:
                    # Shapes of RHS matrices is (N_states, 1), find index with .index(name)
                    state_index, _ = rhs.index(parameter.symbol)
                    T_i = np.take(self.posterior[str(lhs)], state_index, axis=STATE_AXIS)
                    num_i, denom_i = T_i.sum(), T_i.size

                    num += num_i
                    denom += denom_i
            parameters[parameter.name] = num / denom

        return parameters

GMMOptimizer

Bases: EMOptimizer

optimizes parameter values of GMM (sub) model doesnt use guess directly but instead finds next values for parmater from posterior probabilities

rhs of model should all be GMM CompositeNumExprs

Source code in slimfit/minimizers.py

class GMMOptimizer(EMOptimizer):
    """optimizes parameter values of GMM (sub) model
    doesnt use guess directly but instead finds next values for parmater from posterior probabilities

    rhs of model should all be GMM CompositeNumExprs

    """

    # TODO create `step` method which only does one step', execute should be full loop

    def step(self) -> dict[str, float]:
        parameters = {}  # output parameters dictionary

        # All symbols in all 'mu' expressions, then take intersection
        mu_symbols = set.union(*(rhs["mu"].symbols for rhs in self.model.values()))
        mu_symbols &= self.free_parameters.symbols

        # take only the mu symbos in the set of symbols designated as parameters
        mu_parameters = mu_symbols
        # mu_parameters = reduce(
        #     or_, [rhs["mu"].free_parameters.keys() for rhs in self.model.values()]
        # )
        for mu_symbol in mu_parameters:
            num, denom = 0.0, 0.0
            for lhs, gmm_rhs in self.model.items():
                # check if the current mu parameter in this GMM
                if mu_symbol in gmm_rhs["mu"].symbols:
                    col, state_index = gmm_rhs["mu"].index(mu_symbol)
                    T_i = np.take(self.posterior[str(lhs)], state_index, axis=STATE_AXIS)

                    # independent data should be given in the same shape as T_i
                    # which is typically (N, 1), to be sure shapes match we reshape independent data
                    # data = self.xdata[gmm_rhs['x'].name]
                    num += np.sum(T_i * self.xdata[gmm_rhs["x"].name].reshape(T_i.shape))
                    denom += np.sum(T_i)

            parameters[mu_symbol.name] = num / denom

        sigma_symbols = set.union(*(rhs["sigma"].symbols for rhs in self.model.values()))
        sigma_symbols &= self.free_parameters.symbols

        # sigma_parameters = reduce(
        #     or_, [rhs["sigma"].free_parameters.keys() for rhs in self.model.values()]
        # )
        for sigma_symbol in sigma_symbols:
            num, denom = 0.0, 0.0
            # LHS in numerical models are `str` (at the moment)
            for lhs, gmm_rhs in self.model.items():
                # check if the current sigma parameter in this GMM
                if sigma_symbol in gmm_rhs["sigma"].symbols:
                    col, state_index = gmm_rhs["sigma"].index(sigma_symbol)

                    T_i = np.take(self.posterior[str(lhs)], state_index, axis=STATE_AXIS)

                    # Indexing of the MatrixExpr returns elements of its expr
                    mu_name: str = gmm_rhs["mu"][col, state_index].name

                    # Take the corresponding value from the current parameters dict, if its not
                    # there, it must be in the fixed parameters of the model
                    try:
                        mu_value = parameters[mu_name]
                    except KeyError:
                        mu_value = self.fixed_parameters.guess[mu_name]

                    num += np.sum(
                        T_i * (self.xdata[gmm_rhs["x"].name].reshape(T_i.shape) - mu_value) ** 2
                    )

                    denom += np.sum(T_i)

            parameters[sigma_symbol.name] = np.sqrt(num / denom)

        return parameters

LikelihoodOptimizer

Bases: Minimizer

Assumed loss is LogLoss

Source code in slimfit/minimizers.py

class LikelihoodOptimizer(Minimizer):
    """
    Assumed `loss` is `LogLoss`

    """

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        if not isinstance(self.loss, LogSumLoss):
            warnings.warn("Using non-log loss in likelihood optimizer")
        param_shapes = {p.name: p.shape for p in self.free_parameters}
        # TODO rename and generalize?
        self.objective = ScipyObjective(
            model=self.model.numerical,
            loss=self.loss,
            xdata=self.xdata | self.fixed_parameters.guess,
            ydata=self.ydata,
            shapes=param_shapes,
        )

    def execute(
        self, max_iter=250, patience=5, stop_loss=1e-7, verbose=True, **kwargs
    ) -> FitResult:
        """
        kwargs: additional kwargs for scipy minimizer part

        """

        # parameters which needs to be passed / inferred
        # Split top-level multiplications in the model as they can be optimized in log likelihood independently
        # TODO model should have this as property / method
        components: list[tuple[Symbol, NumExprBase]] = []  # todo tuple LHS as variable
        for lhs, rhs in self.model.numerical.items():
            if isinstance(rhs, Mul):
                components += [(lhs, elem) for elem in rhs.values()]
            else:
                components.append((lhs, rhs))

        # Find the sets of components which have common parameters and thus need to be optimized together
        sub_models = intersecting_component_symbols(components, self.free_parameters.symbols)

        # Remove sub models without symbols
        sub_models = [m for m in sub_models if m.symbols]

        pbar = trange(max_iter, disable=not verbose)
        t0 = time.time()

        parameters_current = self.free_parameters.guess  # initialize parameters
        prev_loss = 0.0
        no_progress = 0
        base_result = {}
        for i in pbar:
            y_model = self.model(**parameters_current, **self.xdata, **self.fixed_parameters.guess)
            loss = self.loss(self.ydata, y_model)

            # posterior dict has values with shapes equal to y_model
            # which is (dataopints, states, 1)
            posterior = {k: v / v.sum(axis=STATE_AXIS, keepdims=True) for k, v in y_model.items()}

            # dictionary of new parameter values for in this iteration
            parameters_step = {}
            common_kwargs = dict(
                loss=self.loss,
                xdata=self.xdata,
                ydata=self.ydata,
                posterior=posterior,
            )
            for sub_model in sub_models:
                # At the moment we assume all callables in the sub models to be MatrixCallables
                # determine the kind
                kinds = [getattr(c, "kind", None) for c in sub_model.values()]
                # Filter the general list of parameters to reduce it down to the parameters
                # this model accepts
                # the sub model also needs to the fixed parameters to correctly be able to
                # call the model; GMM also needs it for finding sigma
                sub_parameters = sub_model.filter_parameters(self.parameters)
                if all(k == "constant" for k in kinds):
                    # Constant optimizer doesnt use loss, xdata, ydata,
                    opt = ConstantOptimizer(sub_model, sub_parameters, **common_kwargs)
                    parameters = opt.step()
                elif all(k == "gmm" for k in kinds) and not sub_parameters.has_bounds:
                    # Loss is not used for GMM optimizer step
                    opt = GMMOptimizer(sub_model, sub_parameters, **common_kwargs)
                    parameters = opt.step()
                else:
                    # Previous code:
                    # updated_parameters = [
                    #     Parameter(**(asdict(p) | {"guess": parameters_current.get(p.name) or self.fixed_parameters.guess[p.name]  }))
                    #     for p in sub_parameters
                    # ]

                    # New; needs improvement as it accesses `_names`
                    current_sub = {
                        k: v for k, v in parameters_current.items() if k in sub_parameters._names
                    }
                    updated_parameters = sub_parameters.update_guess(current_sub)

                    # todo loss is not used; should be EM loss while the main loop uses Log likelihood loss
                    opt = ScipyEMOptimizer(sub_model, updated_parameters, **common_kwargs)

                    scipy_result = opt.execute(**kwargs)
                    parameters = scipy_result.fit_parameters
                    base_result["scipy"] = scipy_result

                # collect parameters of this sub_model into parameters dict
                parameters_step |= parameters

            # update for next iteration
            parameters_current = parameters_step

            # loss
            improvement = prev_loss - loss
            prev_loss = loss
            pbar.set_postfix({"improvement": improvement})

            if np.isnan(improvement):
                break
            elif improvement < stop_loss:
                no_progress += 1
            else:
                no_progress = 0

            if no_progress > patience:
                break

        tdelta = time.time() - t0
        gof_qualifiers = {
            "loss": loss,
            "log_likelihood": -loss,
            "n_iter": i + 1,
            "elapsed": tdelta,
            "iter/s": (i + 1) / tdelta,
        }

        result = FitResult(
            fit_parameters=parameters_current,
            fixed_parameters=self.fixed_parameters.guess,
            gof_qualifiers=gof_qualifiers,
            guess=self.free_parameters.guess,
            base_result=base_result,
        )

        return result

execute(max_iter=250, patience=5, stop_loss=1e-07, verbose=True, **kwargs)

kwargs: additional kwargs for scipy minimizer part

Source code in slimfit/minimizers.py

def execute(
    self, max_iter=250, patience=5, stop_loss=1e-7, verbose=True, **kwargs
) -> FitResult:
    """
    kwargs: additional kwargs for scipy minimizer part

    """

    # parameters which needs to be passed / inferred
    # Split top-level multiplications in the model as they can be optimized in log likelihood independently
    # TODO model should have this as property / method
    components: list[tuple[Symbol, NumExprBase]] = []  # todo tuple LHS as variable
    for lhs, rhs in self.model.numerical.items():
        if isinstance(rhs, Mul):
            components += [(lhs, elem) for elem in rhs.values()]
        else:
            components.append((lhs, rhs))

    # Find the sets of components which have common parameters and thus need to be optimized together
    sub_models = intersecting_component_symbols(components, self.free_parameters.symbols)

    # Remove sub models without symbols
    sub_models = [m for m in sub_models if m.symbols]

    pbar = trange(max_iter, disable=not verbose)
    t0 = time.time()

    parameters_current = self.free_parameters.guess  # initialize parameters
    prev_loss = 0.0
    no_progress = 0
    base_result = {}
    for i in pbar:
        y_model = self.model(**parameters_current, **self.xdata, **self.fixed_parameters.guess)
        loss = self.loss(self.ydata, y_model)

        # posterior dict has values with shapes equal to y_model
        # which is (dataopints, states, 1)
        posterior = {k: v / v.sum(axis=STATE_AXIS, keepdims=True) for k, v in y_model.items()}

        # dictionary of new parameter values for in this iteration
        parameters_step = {}
        common_kwargs = dict(
            loss=self.loss,
            xdata=self.xdata,
            ydata=self.ydata,
            posterior=posterior,
        )
        for sub_model in sub_models:
            # At the moment we assume all callables in the sub models to be MatrixCallables
            # determine the kind
            kinds = [getattr(c, "kind", None) for c in sub_model.values()]
            # Filter the general list of parameters to reduce it down to the parameters
            # this model accepts
            # the sub model also needs to the fixed parameters to correctly be able to
            # call the model; GMM also needs it for finding sigma
            sub_parameters = sub_model.filter_parameters(self.parameters)
            if all(k == "constant" for k in kinds):
                # Constant optimizer doesnt use loss, xdata, ydata,
                opt = ConstantOptimizer(sub_model, sub_parameters, **common_kwargs)
                parameters = opt.step()
            elif all(k == "gmm" for k in kinds) and not sub_parameters.has_bounds:
                # Loss is not used for GMM optimizer step
                opt = GMMOptimizer(sub_model, sub_parameters, **common_kwargs)
                parameters = opt.step()
            else:
                # Previous code:
                # updated_parameters = [
                #     Parameter(**(asdict(p) | {"guess": parameters_current.get(p.name) or self.fixed_parameters.guess[p.name]  }))
                #     for p in sub_parameters
                # ]

                # New; needs improvement as it accesses `_names`
                current_sub = {
                    k: v for k, v in parameters_current.items() if k in sub_parameters._names
                }
                updated_parameters = sub_parameters.update_guess(current_sub)

                # todo loss is not used; should be EM loss while the main loop uses Log likelihood loss
                opt = ScipyEMOptimizer(sub_model, updated_parameters, **common_kwargs)

                scipy_result = opt.execute(**kwargs)
                parameters = scipy_result.fit_parameters
                base_result["scipy"] = scipy_result

            # collect parameters of this sub_model into parameters dict
            parameters_step |= parameters

        # update for next iteration
        parameters_current = parameters_step

        # loss
        improvement = prev_loss - loss
        prev_loss = loss
        pbar.set_postfix({"improvement": improvement})

        if np.isnan(improvement):
            break
        elif improvement < stop_loss:
            no_progress += 1
        else:
            no_progress = 0

        if no_progress > patience:
            break

    tdelta = time.time() - t0
    gof_qualifiers = {
        "loss": loss,
        "log_likelihood": -loss,
        "n_iter": i + 1,
        "elapsed": tdelta,
        "iter/s": (i + 1) / tdelta,
    }

    result = FitResult(
        fit_parameters=parameters_current,
        fixed_parameters=self.fixed_parameters.guess,
        gof_qualifiers=gof_qualifiers,
        guess=self.free_parameters.guess,
        base_result=base_result,
    )

    return result

parse_scipy_options(minimizer_options)

rename options for scipy minimizer

Source code in slimfit/minimizers.py

def parse_scipy_options(minimizer_options: dict[str, Any]) -> dict[str, Any]:
    """rename options for scipy minimizer"""

    kwargs = minimizer_options.copy()
    options = minimizer_options.get("options", {})
    if "max_iter" in kwargs:
        options["maxiter"] = kwargs.pop("max_iter")
    if "disp" in kwargs:
        options["disp"] = kwargs.pop("disp")
    if options:
        kwargs["options"] = options

    return kwargs