Example 2 — Custom Python growth model

Demonstrates defining a growth model as a Python function and fitting it via NLModel.

Full script: examples/scripts/03_custom_model.py

Define a Gompertz model

import numpy as np
from pykinbiont import NLModel

def gompertz(p, t):
    """
    Modified Gompertz function.

    Parameters
    ----------
    p : array-like, shape (3,)
        [A, mu, lam]  —  carrying capacity, max growth rate, lag time
    t : float or ndarray
        Time point(s)
    """
    A, mu, lam = p[0], p[1], p[2]
    return A * np.exp(-np.exp(mu * np.e / A * (lam - t) + 1))

custom_gompertz = NLModel(
    name="gompertz",
    func=gompertz,
    param_names=["A", "mu", "lam"],
)

Synthetic data with a visible lag phase

import pykinbiont
from pykinbiont import GrowthData

pykinbiont.configure("/path/to/KinBiont.jl")

def true_gompertz(t, A=1.1, mu=0.5, lam=2.0):
    return A * np.exp(-np.exp(mu * np.e / A * (lam - t) + 1))

rng   = np.random.default_rng(42)
times = np.linspace(0, 18, 60)

curves = np.stack([
    np.maximum(
        true_gompertz(times, A=1.0 + 0.1 * i, mu=0.4 + 0.05 * i, lam=2.0 + 0.3 * i)
        + rng.normal(0, 0.007, len(times)),
        1e-4,
    )
    for i in range(3)
])

data = GrowthData(curves=curves, times=times, labels=["G1", "G2", "G3"])

Fit

from pykinbiont import ModelSpec, FitOptions, fit

spec = ModelSpec(
    models=[custom_gompertz],
    params=[[1.0, 0.5, 2.0]],
    lower=[[0.0, 0.0, 0.0]],
    upper=[[3.0, 3.0, 10.0]],
)
opts = FitOptions(multistart=True, n_restart=30)

results = fit(data, spec, opts)

for r in results:
    p = dict(zip(r.param_names, r.best_params))
    print(f"{r.label}:  A={p['A']:.3f}  mu={p['mu']:.3f}  lam={p['lam']:.3f}"
          f"  AIC={r.best_aic:.2f}")

Compare custom vs built-in

from pykinbiont import MODEL_REGISTRY

logistic_model = MODEL_REGISTRY["NL_logistic"]

spec_both = ModelSpec(
    models=[custom_gompertz, logistic_model],
    params=[[1.0, 0.5, 2.0], [1.2, 0.5, 0.01]],
    lower=[[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]],
    upper=[[3.0, 3.0, 10.0], [5.0, 5.0, 1.0]],
)

results_both = fit(data, spec_both, opts)
df = results_both.to_dataframe()
print(df[["label", "best_model", "aic", "loss"]])

Custom ODE model

For models defined as differential equations, use ODEModel:

from pykinbiont import ODEModel

def logistic_ode(du, u, p, t):
    """Logistic growth ODE: dN/dt = mu*N*(1 - N/K)"""
    mu, K = p[0], p[1]
    du[0] = mu * u[0] * (1.0 - u[0] / K)

ode_model = ODEModel(
    name="logistic_ode",
    func=logistic_ode,
    param_names=["mu", "K"],
    n_eq=1,
)

spec_ode = ModelSpec(
    models=[ode_model],
    params=[[0.5, 1.2]],
    lower=[[0.0, 0.1]],
    upper=[[5.0, 5.0]],
)

results_ode = fit(data, spec_ode, opts)