Examples

Minimizing Python functions

zfit minimizers can also be used to minimize normal Python functions as shown here.

"""This example illustrates how to minimize an arbitrary Python function using the zfit minimizer.

This may has some overhead in the beginning and won't be instantly fast compared to other libraries if run once.

Copyright (c) 2021 zfit
"""

#  Copyright (c) 2024 zfit
from __future__ import annotations

import numpy as np

import zfit

# set everything to numpy mode. This is only needed if we
# don't use
zfit.run.set_autograd_mode(False)
zfit.run.set_graph_mode(False)

# create our favorite minimizer
minimizer = zfit.minimize.IpyoptV1()


# minimizer = zfit.minimize.Minuit()
# minimizer = zfit.minimize.ScipyTrustConstrV1()
# minimizer = zfit.minimize.NLoptLBFGSV1()


def func(x):
    x = np.asarray(x)  # make sure it's an array
    return np.sum((x - 0.1) ** 2 + x[1] ** 4)


# We can use the same in pure TF, then we can also use
# the analytic gradient
#
# import tensorflow as tf
#
# @tf.function
# def func(x):
#     x = tf.convert_to_tensor(x)  # make sure it's an array
#     return tf.reduce_sum((x - 0.1) ** 2 + x[1] ** 4)


# we can also use a more complicated function instead
# from scipy.optimize import rosen as func


# we need to set the errordef, the definition of "1 sigma"
func.errordef = 0.5

# initial parameters
params = [1, -3, 2, 1.4, 11]
# or for a more fine-grained control
# params = {
#     'value': [1, -3, 2, 1.4, 11],  # mandatory
#     'lower': [-2, -5, -5, -10, -15],  # lower bound, can be omitted
#     'upper': [2, 4, 5, 10, 15],  # upper bound, can be omitted
#     'step_size': [0.1] * 5,  # initial step size, can be omitted
# }

# minimize
result = minimizer.minimize(func, params)

# estimate errors
result.hesse()
result.errors()

Composing PDFs

#  Copyright (c) 2024 zfit
from __future__ import annotations

import zfit

# create space
obs = zfit.Space("x", -10, 10)

# parameters
mu = zfit.Parameter("mu", 1.0, -4, 6)
sigma = zfit.Parameter("sigma", 1.0, 0.1, 10)
lambd = zfit.Parameter("lambda", -1.0, -5.0, 0)
frac = zfit.Parameter("fraction", 0.5, 0.0, 1.0)

# pdf creation
gauss = zfit.pdf.Gauss(mu=mu, sigma=sigma, obs=obs)
exponential = zfit.pdf.Exponential(lambd, obs=obs)

sum_pdf = zfit.pdf.SumPDF([gauss, exponential], fracs=frac)

free_param = zfit.Parameter("free", 5, 0, 15)
param_shift = zfit.ComposedParameter("comp", lambda x: x + 5, params=free_param)

Custom 3D Functor

#  Copyright (c) 2024 zfit

# Functors are PDfs and Functions that depend on other PDFs or Functions. They can be used to define
# in a custom way combinations of PDFs or wrapping a single PDF.
# An example would be to create the sum of two PDFs. Of course this is already implemented in zfit
# as the functor SumPDF([pdf1, pdf2], fracs=...). For advanced uses, you
# can define your own functor as demonstrated below.
# ruff: noqa: F821
from __future__ import annotations

import numpy as np

import zfit
import zfit.z


class CombinePolynomials(zfit.pdf.BaseFunctor):
    """Example of a functor pdf that adds three polynomials.

    DEMONSTRATION PURPOSE ONLY, DO **NOT** USE IN REAL CASE.
    """

    def __init__(self, angle1, angle2, angle3, name="3dPolynomial"):
        pdfs = [angle1, angle2, angle3]
        super().__init__(pdfs=pdfs, name=name)

    _N_OBS = 3

    def _unnormalized_pdf(self, x):
        cosangle1_data, cosangle2_data, angle3_data = z.unstack_x(x)
        cosangle1, cosangle2, angle3 = self.pdfs

        return cosangle1.pdf(cosangle1_data) + cosangle2.pdf(cosangle2_data) + angle3.pdf(angle3_data)


def create_angular():
    # create parameters
    # c00 = zfit.Parameter(...)
    # ... and so on

    cosangle1_space = zfit.Space("cos angle 1", limits=(-1, 1))
    cosangle2_space = zfit.Space("cos angle 2", limits=(-1, 1))
    angle3_space = zfit.Space("angle 3", limits=(-np.pi, np.pi))

    # this part could also be moved inside the __init__, but then the init would need to take all
    # the coefficients and the spaces as arguments
    cosangle1_pdf = zfit.pdf.Chebyshev(obs=cosangle1_space, coeffs=[c00, c01, c02])
    cosangle2_pdf = zfit.pdf.Chebyshev(obs=cosangle2_space, coeffs=[c10, c11, c12])
    angle3_pdf = zfit.pdf.Chebyshev(obs=angle3_space, coeffs=[c20, c21, c22])

    return CombinePolynomials(angle1=cosangle1_pdf, angle2=cosangle2_pdf, angle3=angle3_pdf)

Custom PDF simple

#  Copyright (c) 2024 zfit
from __future__ import annotations

import zfit
import zfit.z.numpy as znp


class CustomPDF(zfit.pdf.ZPDF):
    """1-dimensional PDF implementing the exp(alpha * x) shape."""

    _PARAMS = ("alpha",)  # specify which parameters to take

    def _unnormalized_pdf(self, x):  # implement function
        data = x[0]  # axis 0
        alpha = self.params["alpha"]

        return znp.exp(alpha * data)


obs = zfit.Space("obs1", -4, 4)

custom_pdf = CustomPDF(obs=obs, alpha=0.2)

integral = custom_pdf.integrate(limits=(-1, 2))
sample = custom_pdf.sample(n=1000)
prob = custom_pdf.pdf(sample)

Custom PDF advanced

#  Copyright (c) 2024 zfit
from __future__ import annotations

import tensorflow as tf

import zfit
from zfit import z


class CustomPDF2D(zfit.pdf.BasePDF):
    """My custom, 2 dimensional pdf where the axes are: Energy, Momentum."""

    def __init__(
        self,
        param1,
        param2,
        param3,
        obs,
        name="CustomPDF",
    ):
        # we can now do complicated stuff here if needed
        # only thing: we have to specify explicitly here what is which parameter
        params = {
            "super_param": param1,  # we can change/compose etc parameters
            "param2": param2,
            "param3": param3,
        }
        super().__init__(obs, params, name=name)

    @zfit.supports()
    def _unnormalized_pdf(self, x, params):
        energy = x[0]
        momentum = x[1]
        param1 = params["super_param"]
        param2 = params["param2"]
        param3 = params["param3"]

        # just a fantasy function
        return param1 * tf.cos(energy**2) + tf.math.log(param2 * momentum**2) + param3


# add an analytic integral


# define the integral function
def integral_full(limits, norm_range, params, model):
    del norm_range, model  # not used here
    lower, upper = limits.v1.limits
    param1 = params["super_param"]
    param2 = params["param2"]
    param3 = params["param3"]

    # calculate the integral here, dummy integral, wrong!
    return param1 * param2 * param3 + z.reduce_sum([lower, upper])


# define the space over which it is defined. Here, we use the axes
integlimits_axis0 = zfit.Space(axes=0, lower=-10, upper=10)
integlimits_axis1 = zfit.Space(axes=1, lower=zfit.Space.ANY_LOWER, upper=zfit.Space.ANY_UPPER)

integral_full_limits = integlimits_axis0 * integlimits_axis1  # creates 2D space
CustomPDF2D.register_analytic_integral(func=integral_full, limits=integral_full_limits)


# define the partial integral function
def integral_axis1(x, limits, norm_range, params, model):
    del norm_range, model  # not used here

    data_0 = x[0]  # data from axis 0

    param1 = params["super_param"]
    param2 = params["param2"]
    param3 = params["param3"]

    lower, upper = limits.limit1d  # for a more detailed guide, see the space.py example
    lower = z.convert_to_tensor(lower)  # the limits are now 1-D, for axis 1
    upper = z.convert_to_tensor(upper)

    # calculate the integral here, dummy integral
    return data_0**2 * param1 * param2 * param3 + z.reduce_sum([lower, upper])
    # notice that the returned shape will be in the same as data_0, e.g. the number of events given in x


# define the space over which it is defined. Here, we use the axes
integral_axis1_limits = zfit.Space(
    axes=(1,),  # axes one corresponds to the second obs, here obs2
    lower=zfit.Space.ANY_LOWER,
    upper=zfit.Space.ANY_UPPER,
)

CustomPDF2D.register_analytic_integral(func=integral_axis1, limits=integral_axis1_limits)

if __name__ == "__main__":
    import numpy as np

    obs = zfit.Space("obs1", -10, 10) * zfit.Space("obs2", -3, 5)
    pdf = CustomPDF2D(
        param1=1, param2=2, param3=3, obs=obs
    )  # if a Python number is passed, it will be regarded as a constant
    sample = pdf.sample(n=1000)
    pdf.pdf([[2.0, 2.5], [5.4, 3.2]])
    # x_part = zfit.Data(np.array([2.1, 2.2, 3.2]), obs="obs1")
    x_part = np.array([2.1, 2.2, 3.2])
    # integrate over obs2 with limits 1, 2 for the `x_part`. This will use the analytic integral above
    pdf.partial_integrate(x=x_part, limits=zfit.Space("obs2", 1, 2))
    # we can explicitly call the analytic integral. Without registering it (e.g. comment the line with the `register`
    # and run again), it will raise an error
    pdf.partial_analytic_integrate(x=x_part, limits=zfit.Space("obs2", 1, 2))

Custom 3D PDF

#  Copyright (c) 2024 zfit
from __future__ import annotations

import tensorflow as tf

import zfit

# IMPORTANT! The communication of which axis corresponds to which data point happens here. So the user knows now that
# he should create this pdf with a space in the obs (x, y, z).


class CustomPDF(zfit.pdf.ZPDF):
    """3-dimensional PDF calculating doing something fancy, takes x, y, z as data."""

    _PARAMS = ("alpha", "beta")  # specify which parameters to take
    _N_OBS = 3

    @zfit.supports()
    def _unnormalized_pdf(self, x, params):  # implement function
        x0 = x[0]
        x1 = x[1]
        x2 = x[2]
        alpha = params["alpha"]
        beta = params["beta"]
        x0_new = tf.math.cos(alpha) * x0
        x1_new = tf.math.sinh(beta) * x1
        x2_new = x2 + 4.2
        return x0_new**2 + x1_new**2 + x2_new**2


xobs = zfit.Space("xobs", -4, 4)
yobs = zfit.Space("yobs", -3, 3)
zobs = zfit.Space("z", -2, 2)
obs = xobs * yobs * zobs

alpha = zfit.Parameter("alpha", 0.2)  # floating
beta = zfit.Parameter("beta", 0.4, floating=False)  # non-floating
custom_pdf = CustomPDF(obs=obs, alpha=alpha, beta=beta)

integral = custom_pdf.integrate(limits=obs)  # = 1 since normalized
sample = custom_pdf.sample(n=1000)  # DO NOT USE THIS FOR TOYS!
prob = custom_pdf.pdf(sample)  # DO NOT USE THIS FOR TOYS!

Multidimensional preprocess fit

#  Copyright (c) 2024 zfit
from __future__ import annotations

import numpy as np

import zfit

# create space
xobs = zfit.Space("xobs", -4, 4)
yobs = zfit.Space("yobs", -3, 5)
zobs = zfit.Space("z", -2, 4)
obs = xobs * yobs * zobs

# parameters
mu1 = zfit.Parameter("mu1", 1.0, -4, 6)
mu23 = zfit.Parameter("mu_shared", 1.0, -4, 6)
sigma12 = zfit.Parameter("sigma_shared", 1.0, 0.1, 10)
sigma3 = zfit.Parameter("sigma3", 1.0, 0.1, 10)

# model building, pdf creation
gauss_x = zfit.pdf.Gauss(mu=mu1, sigma=sigma12, obs=xobs)
gauss_y = zfit.pdf.Gauss(mu=mu23, sigma=sigma12, obs=yobs)
gauss_z = zfit.pdf.Gauss(mu=mu23, sigma=sigma3, obs=zobs)

product_gauss = zfit.pdf.ProductPDF([gauss_x, gauss_y, gauss_z])

# data
normal_np = np.random.normal(loc=[2.0, 2.5, 2.5], scale=[3.0, 3, 1.5], size=(10000, 3))
data_raw = zfit.Data(normal_np, obs=obs)  # or from anywhere else, e.g. root

df = data_raw.to_pandas()
# preprocessing here, rename things. Match column names with the observable names "xobs", "yobs", "z" (they have to be
# contained, more columns in the df is not a problem)

data = df  # we can directly use the dataframe
# data = normal_np  # or the numpy array
# data = zfit.Data.from_pandas(df, obs=obs)  # or create another zfit data object

# create NLL
nll = zfit.loss.UnbinnedNLL(model=product_gauss, data=data)

# create a minimizer
minimizer = zfit.minimize.Minuit()
result = minimizer.minimize(nll).update_params()

# do the error calculations, here with minos
param_errors, _ = result.errors()

Plot signal background

#  Copyright (c) 2024 zfit
from __future__ import annotations

import matplotlib.pyplot as plt
import mplhep
import numpy as np

import zfit

mplhep.style.use("LHCb2")
# create space
obs = zfit.Space("x", -10, 10)

# parameters
mu = zfit.Parameter("mu", 1.0, -4, 6)
sigma = zfit.Parameter("sigma", 1.0, 0.1, 10)
lambd = zfit.Parameter("lambda", -0.06, -1, -0.01)
frac = zfit.Parameter("fraction", 0.3, 0, 1)

# model building, pdf creation
gauss = zfit.pdf.Gauss(mu=mu, sigma=sigma, obs=obs)
exponential = zfit.pdf.Exponential(lambd, obs=obs)
model = zfit.pdf.SumPDF([gauss, exponential], fracs=frac)

# for a quick plot, we can use the attached "plotter"
plt.figure()
plt.title("Model before fit, using plotter")
model.plot.plotpdf()

# data
n_sample = 5000

exp_data = exponential.sample(n=n_sample * (1 - frac))

gauss_data = gauss.sample(n=n_sample * frac)

data = model.create_sampler(n_sample, limits=obs)
data.resample()

# set the values of the parameters
zfit.param.set_values([mu, sigma, lambd, frac], [0.5, 1.2, -0.05, 0.07])

# alternatively, we can set the values individually
# mu.set_value(0.5)
# sigma.set_value(1.2)
# lambd.set_value(-0.05)
# frac.set_value(0.07)

# plot the data
n_bins = 50

plot_scaling = n_sample / n_bins * obs.volume

x = np.linspace(obs.v1.lower, obs.v1.upper, 1000)


def plot_pdf(title):
    plt.figure()
    plt.title(title)
    y = model.pdf(x)
    y_gauss = gauss.pdf(x) * frac
    y_exp = exponential.pdf(x) * (1 - frac)
    plt.plot(x, y * plot_scaling, label="Sum - Model")
    plt.plot(x, y_gauss * plot_scaling, label="Gauss - Signal")
    plt.plot(x, y_exp * plot_scaling, label="Exp - Background")
    mplhep.histplot(
        data.to_binned(n_bins),
        yerr=True,
        color="black",
        histtype="errorbar",
    )
    plt.ylabel("Counts")
    plt.xlabel("$B_{mass}$")


# plot the pdf BEFORE fitting
plot_pdf("Before fitting")


# create NLL
nll = zfit.loss.UnbinnedNLL(model=model, data=data)

# create a minimizer
minimizer = zfit.minimize.Minuit()
result = minimizer.minimize(nll).update_params()

# do the error calculations, here with minos
param_errors, _ = result.errors()

plot_pdf(title="After fitting")
plt.legend()

# uncomment to display plots
# plt.show()

Simple fit

#  Copyright (c) 2024 zfit
from __future__ import annotations

import numpy as np

import zfit

# create space
obs = zfit.Space("x", -2, 3)

# parameters
mu = zfit.Parameter("mu", 1.2, -4, 6)
sigma = zfit.Parameter("sigma", 1.3, 0.5, 10)

# model building, pdf creation
gauss = zfit.pdf.Gauss(mu=mu, sigma=sigma, obs=obs)

# data
data = np.random.normal(loc=2.0, scale=3.0, size=10000)
# we can convert it to a zfit Data object to make sure the data is within the limits or use the space to filter manually
data = obs.filter(data)  # works also for pandas DataFrame
# data = zfit.Data.from_numpy(obs=obs, array=data)

# create NLL
nll = zfit.loss.UnbinnedNLL(model=gauss, data=data)

# create a minimizer
minimizer = zfit.minimize.Minuit(gradient=False)
result = minimizer.minimize(nll).update_params()

# do the error calculations with a hessian approximation
param_errors = result.hesse()

# or here with minos
param_errors_asymmetric, new_result = result.errors()

Signal-background-mass fit

#  Copyright (c) 2024 zfit
from __future__ import annotations

import zfit

# create space
obs = zfit.Space("x", -10, 10)

# parameters
mu = zfit.Parameter("mu", 1.0, -4, 6)
sigma = zfit.Parameter("sigma", 1.0, 0.1, 10)
lambd = zfit.Parameter("lambda", -0.06, -1, -0.01)
frac = zfit.Parameter("fraction", 0.3, 0, 1)

# model building, pdf creation
gauss = zfit.pdf.Gauss(mu=mu, sigma=sigma, obs=obs)
exponential = zfit.pdf.Exponential(lambd, obs=obs)
model = zfit.pdf.SumPDF([gauss, exponential], fracs=frac)

# data
n_sample = 10000

# if we sample once, otherwise use create_sampler
data = model.sample(n_sample, limits=obs)  # limits can be omitted, then the default limits are used

# set the values to a start value for the fit
zfit.param.set_values({mu: 0.5, sigma: 1.2, lambd: -0.05, frac: 0.07})

# create NLL
nll = zfit.loss.UnbinnedNLL(model=model, data=data)

# create a minimizer
minimizer = zfit.minimize.Minuit()
zfit.run.experimental_disable_param_update()
result = minimizer.minimize(
    nll
).update_params()  # forward compatibility, update_params will update the values of the parameters to the minimum

# do the error calculations, here with minos
param_hesse = result.hesse()
param_errors, new_result = result.errors()

Signal-background-mass fit extended

#  Copyright (c) 2024 zfit
from __future__ import annotations

import pickle

import zfit

n_bins = 50

# create space
obs = zfit.Space("x", -10, 10)

# parameters
mu = zfit.Parameter("mu", 1.0, -4, 6)
sigma = zfit.Parameter("sigma", 1.0, 0.1, 10)
lambd = zfit.Parameter("lambda", -0.06, -1, -0.01)
n_bkg = zfit.Parameter("n_bkg", 20000)
n_sig = zfit.Parameter("n_sig", 1000)

# model building, pdf creation
gauss_extended = zfit.pdf.Gauss(mu=mu, sigma=sigma, obs=obs, extended=n_sig)
exp_extended = zfit.pdf.Exponential(lambd, obs=obs, extended=n_bkg)

model = zfit.pdf.SumPDF([gauss_extended, exp_extended])

data = model.create_sampler(n=21200)
data.resample()

# set the values to a start value for the fit
zfit.param.set_values({mu: 0.5, sigma: 1.2, lambd: -0.05})
# alternatively, we can set the values individually
# mu.set_value(0.5)
# sigma.set_value(1.2)
# lambd.set_value(-0.05)

# create NLL
nll = zfit.loss.ExtendedUnbinnedNLL(model=model, data=data)

# create a minimizer
# this uses the automatic gradient (more precise, but can fail)
minimizer = zfit.minimize.Minuit(gradient="zfit")
result = minimizer.minimize(nll).update_params()
# do the error calculations, here with hesse, than with minos
param_hesse = result.hesse()

# the second return value is a new FitResult if a new minimum was found

(
    param_errors,
    _,
) = result.errors()

# Storing the result can be achieved in many ways. Using dill (like pickle), we can dump and load the result
result_dilled = zfit.dill.dumps(result)
result_loaded = zfit.dill.loads(result_dilled)

# EXPERIMENTAL: we can serialize the model to a human-readable format with HS3
# human readable representation
hs3like = zfit.hs3.dumps(nll)
# print(hs3like)
# and we can load it again
nll_loaded = zfit.hs3.loads(hs3like)

# pickle can also be used, but as not all objects are pickleable, we need to freeze the result
# which makes it read-only and no new error calculations can be done
result.freeze()
dumped = pickle.dumps(result)
loaded = pickle.loads(dumped)

zfit.param.set_values(model.get_params(), loaded)

Simultanous fit

#  Copyright (c) 2024 zfit

# create space
from __future__ import annotations

import numpy as np

import zfit

obs = zfit.Space("x", -10, 10)

# parameters
mu_shared = zfit.Parameter("mu_shared", 1.0, -4, 6)  # mu is a shared parameter
sigma1 = zfit.Parameter("sigma_one", 1.0, 0.1, 10)
sigma2 = zfit.Parameter("sigma_two", 1.0, 0.1, 10)

# model building, pdf creation
gauss1 = zfit.pdf.Gauss(mu=mu_shared, sigma=sigma1, obs=obs)
gauss2 = zfit.pdf.Gauss(mu=mu_shared, sigma=sigma2, obs=obs)

# data
normal_np = np.random.normal(loc=2.0, scale=3.0, size=10000)
data1 = zfit.Data.from_numpy(obs=obs, array=normal_np)  # data for the first model, can also be just a numpy array
# the data objects just makes sure that the data is within the limits

# data
normal_np = np.random.normal(loc=2.0, scale=4.0, size=10000)
data2 = zfit.Data.from_numpy(obs=obs, array=normal_np)

# create simultaenous loss, two possibilities
nll_simultaneous = zfit.loss.UnbinnedNLL(model=[gauss1, gauss2], data=[data1, data2])
# OR, equivalently
nll1 = zfit.loss.UnbinnedNLL(model=gauss1, data=data1)
nll2 = zfit.loss.UnbinnedNLL(model=gauss2, data=data2)
nll_simultaneous2 = nll1 + nll2

minimizer = zfit.minimize.Minuit()
minimizer.minimize(nll_simultaneous2).update_params()

Spaces

#  Copyright (c) 2024 zfit
from __future__ import annotations

import numpy as np

import zfit

# create a multidimensional space
xobs = zfit.Space("xobs", -4, 4)
yobs = zfit.Space("yobs", -3, 5)
zobs = zfit.Space("z", -2, 4)
obs = xobs * yobs * zobs
assert obs.n_obs == 3

# retrieve a subspace
xobs = obs.with_obs("xobs")
yxobs = obs.with_obs(["yobs", "xobs"])  # note the change in order

# retrieve the limits
# they have been overcomplicated in the past, now they are simple. Just access the limits via space.v1
limits = obs.v1.limits  # (lower, upper)
lower = obs.v1.lower
upper = obs.v1.upper

assert np.all(lower == [-4, -3, -2])
assert np.all(upper == [4, 5, 4])

# the volume is the product of the volumes of the individual spaces
volume = obs.volume
assert volume == xobs.volume * yobs.volume * zobs.volume
assert xobs.volume == xobs.v1.upper - xobs.v1.lower

# for example, creating a linspace object is simple
x = np.linspace(xobs.v1.lower, xobs.v1.upper, 1000)
# or
x = np.linspace(*xobs.v1.limits, 1000)

# or even in 3D
x = np.linspace(*obs.v1.limits, 1000)
assert x.shape == (1000, 3)