cgnm 0.9.3.post0

Cluster Gauss-Newton Method – find multiple minimisers of a nonlinear least-squares problem

CGNM – Python Port

Python translation of the Cluster Gauss-Newton Method R package (v0.9.3) by Yasunori Aoki.

Reference: Aoki et al. (2020). Cluster Gauss-Newton method. Optimization and Engineering. https://doi.org/10.1007/s11081-020-09571-2 Aoki and Sugiyama (2023). Cluster Gauss-Newton method for a quick approximation of profile likelihood: With application to physiologically-based pharmacokinetic models. https://doi.org/10.1002/psp4.13055

---

What it does

CGNM finds multiple minimisers of the nonlinear least-squares problem:

$$\operatorname{argmin}_x | f(x) - y^* |^2$$

where f is your model (e.g. a PBPK ODE), y* is observed data, and x are parameters.
It is especially useful when the problem has non-unique solutions (e.g. flip-flop kinetics).

---

Installation

cd python_cgnm
pip install -e .

Dependencies: numpy, scipy, pandas, scikit-learn, matplotlib

Install Streamlit app support:

pip install -e ".[app]"

Launch the app:

cgnm-streamlit

---

Quick start

import numpy as np
import cgnm

observation_time = np.array([0.1, 0.2, 0.4, 0.6, 1, 2, 3, 6, 12])

def model(x):
    ka, V1, CL = x
    t = observation_time
    Cp = ka * 1000 / (V1 * (ka - CL/V1)) * (np.exp(-CL/V1 * t) - np.exp(-ka * t))
    return np.log10(Cp)

observation = np.log10([4.91, 8.65, 12.4, 18.7, 24.3, 24.5, 18.4, 4.66, 0.238])

result = cgnm.cluster_gauss_newton_method(
    nonlinear_function=model,
    target_vector=observation,
    initial_lower_range=[0.1, 0.1, 0.1],
    initial_upper_range=[10.0, 10.0, 10.0],
    num_minimizers_to_find=250,
    num_iteration=25,
    parameter_names=["Ka", "V1", "CL"],
    save_log=False,
)

print(cgnm.accepted_approximate_minimizers(result))

---

Key functions

Function	Description
cluster_gauss_newton_method()	Main optimisation (returns result dict)
cluster_gauss_newton_bootstrap_method()	Bootstrap uncertainty quantification
accepted_approximate_minimizers()	DataFrame of accepted minimisers
accepted_indices()	Integer indices of accepted minimisers
top_indices()	Indices of best-N minimisers
table_parameter_summary()	Quantile summary table
plot_rank_ssr()	SSR vs rank plot
plot_para_distribution_by_histogram()	Parameter distribution histograms
plot_goodness_of_fit()	Goodness-of-fit diagnostics
plot_profile_likelihood()	Profile likelihood (requires save_log=True)
plot_simulation_with_ci()	Model simulation with CI band

---

Differences from the R version

Feature	R	Python
Parallel eval	parallel package (optional)	Pass a vectorised function accepting 2-D input
Log files	.RDATA	pickle (.pkl) in CGNM_log/
Shiny GUI	shinyCGNM()	cgnm-streamlit
Reparameterisation strings	eval(parse(...))	Handled via lower_bound/upper_bound arrays

---

Running the example

python example.py

Produces rank_ssr.png, param_histogram.png, goodness_of_fit.png.

---

ODE Solver Tools

Two companion modules let you define pharmacokinetic / pharmacodynamic models as ODE systems and simulate them — bridging the gap between rxODE (R) and Python.

Module	Purpose
ode_solver.py	Multi-dose ODE integration engine (pure Python / scipy)
rxode_compat.py	Translation helpers — parse rxODE syntax, mirror event-table API

---

ode_solver.py — MultiDoseODESolver

Concept

scipy.integrate.solve_ivp handles one continuous interval at a time. MultiDoseODESolver automatically splits the time axis at every dose event, applies the transition (bolus or infusion), and chains the segments together.

Constructor

from ode_solver import MultiDoseODESolver

solver = MultiDoseODESolver(
    rhs,                          # callable: dxdt = rhs(t, x, params)
    n_states,                     # int: number of state variables
    ic,                           # array-like: initial conditions (all zeros typical)
    dose_times,                   # list of floats: when each dose occurs
    dose_amounts,                 # list of floats: how much drug per dose
    dose_compartments = 0,        # int or list: which state index receives the drug
    dose_durations    = 0.0,      # 0 = bolus; >0 = zero-order infusion of that duration
    set_events        = None,     # list of (time, compartment, value): force state[cpt] = val
    method            = "RK45",   # "RK45" (non-stiff) or "Radau" (stiff)
    rtol              = 1e-6,
    atol              = 1e-9,
)

.solve(obs_times, params, ic=None)

sol = solver.solve(
    obs_times = [1, 4, 12, 24],   # observation times (any order)
    params    = (Ka, V, CL),      # passed as 3rd arg to rhs
)
# sol.t  – observation times, shape (n_obs,)
# sol.y  – state matrix, shape (n_states, n_obs)
# sol.y[state_index, time_index]

Important: amount vs concentration states

If your RHS models state as concentration (e.g. d/dt(C) = -ke*C), pass dose_amounts = Dose / V so the solver adds the correct amount to the concentration state. If your RHS models amount (e.g. d/dt(A) = -ke*A), pass dose_amounts = Dose directly.

RHS state type	dose_amounts to pass
Amount (A)	Dose
Concentration (C = A/V)	Dose / V

Dosing types

# Oral bolus every 24 h for 3 days (depot = state 0)
solver = MultiDoseODESolver(rhs, 2, [0,0],
    dose_times=[0, 24, 48], dose_amounts=[100, 100, 100],
    dose_compartments=0, dose_durations=0)

# 2-hour IV infusion (total 200 mg) into central (state 1)
solver = MultiDoseODESolver(rhs, 2, [0,0],
    dose_times=[0], dose_amounts=[200],
    dose_compartments=1, dose_durations=2.0)

# Reset a compartment to 0 at t=12 (equivalent to rxODE evid=2)
solver = MultiDoseODESolver(rhs, 1, [0],
    dose_times=[0], dose_amounts=[100],
    set_events=[(12.0, 0, 0.0)])

Stiff systems

For stiff ODE systems (e.g. fast/slow compartments, binding kinetics) use method="Radau":

solver = MultiDoseODESolver(rhs, n_states, ic, ..., method="Radau")

---

rxode_compat.py — translation helpers

Purpose: make it easy to migrate existing rxODE / rxode2 R code to Python with minimal changes.

parse_rxode(ode_text) — translate an ODE string

from rxode_compat import parse_rxode, show_rhs

model = parse_rxode("""
    d/dt(A_gut) = -Ka*A_gut
    d/dt(C)     =  Ka*A_gut/V - (CL/V)*C
""")

Equivalent R:

library(rxode2)
RxODE("
    d/dt(A_gut) = -Ka*A_gut
    d/dt(C)     =  Ka*A_gut/V - (CL/V)*C
")

parse_rxode returns a ParsedModel with:

Attribute	Content
.state_names	['A_gut', 'C']
.param_names	['CL', 'Ka', 'V'] (auto-inferred, sorted)
.rhs	Python function rhs(t, x, params) ready for MultiDoseODESolver

Automatic R → Python conversions

R	Python
^	**
exp(x), log(x), sqrt(x)	np.exp(x), np.log(x), np.sqrt(x)
sin, cos, tan, abs	np.sin, np.cos, np.tan, abs
ceiling, floor	np.ceil, np.floor
TRUE / FALSE	True / False
pi, Inf, NaN, NA	np.pi, np.inf, np.nan, np.nan

Algebraic intermediate variables (computed before the ODEs) are supported:

model = parse_rxode("""
    Cp = C / V1          # intermediate — available in ODE expressions below
    d/dt(C) = -Cp * CL
""")

Print the generated Python function for debugging:

show_rhs(model)
# def _generated_rhs(t, x, params):
#     ...

Explicit parameter list (use when auto-detection is ambiguous):

model = parse_rxode(ode_text, param_names=["Ka", "V", "CL"])

---

EventTable and et() — dosing + sampling schedule

Mirrors the rxode2 eventTable() / et() API.

from rxode_compat import et, EventTable

# Single oral dose
ev = et(amt=1000, cmt=1).add_sampling([0.5, 1, 2, 4, 8, 12, 24])

# Multiple oral doses every 12 h
ev = (EventTable()
      .add_dosing(dose=500, nbr_doses=5, dosing_interval=12, dosing_to=1)
      .add_sampling(np.arange(0, 61, 4)))

# 30-min IV infusion (total 200 mg, rate = 400 mg/h) into central (cmt 2)
ev = et().add_dosing(dose=200, dosing_to=2, rate=400).add_sampling([1, 2, 4, 8])

add_dosing parameters

Parameter	R equivalent	Notes
dose	dose=	Amount per event
nbr_doses	nbr.doses=	Number of doses (default 1)
dosing_interval	ii=	Time between doses; required when nbr_doses > 1
start_time	start.time=	Time of first dose (default 0)
dosing_to	dosing.to=	Compartment number 1-indexed (rxODE convention)
rate	rate=	Infusion rate; sets duration = dose / rate
duration	(derived)	Explicit infusion duration; overrides rate
evid	evid=	1 = dose (default), 2 = reset+dose, 3 = reset only

Full rxODE → Python mapping

R (rxode2)                                  Python
──────────────────────────────────────────  ──────────────────────────────────────────
et(amt=100, cmt=1)                          et(amt=100, cmt=1)
add.dosing(dose, nbr.doses=5, ii=12,        .add_dosing(dose, nbr_doses=5,
           dosing.to=1)                                 dosing_interval=12, dosing_to=1)
add.dosing(dose=200, dosing.to=2, rate=400) .add_dosing(dose=200, dosing_to=2, rate=400)
add.sampling(c(0.5, 1, 2, 4))              .add_sampling([0.5, 1, 2, 4])
add.sampling(seq(0, 24, by=0.5))           .add_sampling(np.arange(0, 24.5, 0.5))
evid=2  (reset + dose)                      evid=2
evid=3  (reset only)                        evid=3

---

rxode_solve() — one-call solve → DataFrame

from rxode_compat import parse_rxode, et, rxode_solve

model = parse_rxode("""
    d/dt(A_gut) = -Ka*A_gut
    d/dt(C)     =  Ka*A_gut/V - (CL/V)*C
""")

ev = et(amt=1000, cmt=1).add_sampling([0.5, 1, 2, 4, 8, 12, 24])

result = rxode_solve(model, {"Ka": 1.2, "V": 20, "CL": 3}, ev)
print(result)
#    time        A_gut          C
#     0.5  5.19e+02   21.65
#     1.0  2.85e+02   31.97
#     ...

Equivalent R:

rxSolve(model, list(Ka=1.2, V=20, CL=3), ev)

rxode_solve options:

Parameter	Default	Notes
model	—	ParsedModel from parse_rxode
params	—	dict (preferred) or sequence aligned with model.param_names
events	—	EventTable
method	"RK45"	"Radau" for stiff systems
rtol / atol	1e-6 / 1e-9	Solver tolerances
ic	all zeros	Override initial conditions

---

Complete worked example

One-compartment oral model, 3 doses, compare with R

Python:

from rxode_compat import parse_rxode, et, rxode_solve
import numpy as np

model = parse_rxode("""
    d/dt(A_gut) = -Ka*A_gut
    d/dt(C)     =  Ka*A_gut/V - (CL/V)*C
""")

ev = (et()
      .add_dosing(dose=500, nbr_doses=3, dosing_interval=12, dosing_to=1)
      .add_sampling(np.arange(0, 37, 2)))

result = rxode_solve(model, {"Ka": 0.8, "V": 15, "CL": 2}, ev)
print(result[["time", "C"]])

Equivalent R:

library(rxode2)

mod <- RxODE("
    d/dt(A_gut) = -Ka*A_gut
    d/dt(C)     =  Ka*A_gut/V - (CL/V)*C
")

ev <- et() |>
    add.dosing(dose=500, nbr.doses=3, dosing.interval=12, dosing.to=1) |>
    add.sampling(seq(0, 36, by=2))

rxSolve(mod, list(Ka=0.8, V=15, CL=2), ev)[, c("time", "C")]

---

Using ODE models with CGNM

To use your ODE model as the nonlinear_function for CGNM parameter estimation:

import numpy as np
import cgnm
from rxode_compat import parse_rxode, et, rxode_solve

obs_times   = np.array([0.5, 1, 2, 4, 8, 12, 24])
observation = np.array([...])   # your measured concentrations

model = parse_rxode("""
    d/dt(A_gut) = -Ka*A_gut
    d/dt(C)     =  Ka*A_gut/V - (CL/V)*C
""")

def wrapped_model(x):
    Ka, V, CL = x
    ev = et(amt=1000, cmt=1).add_sampling(obs_times)
    df = rxode_solve(model, {"Ka": Ka, "V": V, "CL": CL}, ev)
    return df["C"].values

result = cgnm.cluster_gauss_newton_method(
    nonlinear_function=wrapped_model,
    target_vector=observation,
    initial_lower_range=[0.1, 1, 0.1],
    initial_upper_range=[5.0, 50, 10.0],
    num_minimizers_to_find=100,
    parameter_names=["Ka", "V", "CL"],
)
print(cgnm.accepted_approximate_minimizers(result))

---

Validation

Run the full numerical validation suite (all results compared against closed-form analytical solutions):

python validate_ode_solver.py

Expected output: 8/8 tests passed, all relative errors < 1e-4.