CardamomOT.inference.degradations¶

Utilities for degradation rate inference and temporal epsilon estimation.

This module provides PyTorch models and helper functions used by the CARDAMOM pipeline when learning gene-specific degradation parameters from protein dynamics. It includes the GeneRegulatoryODE_softmax neural ODE model and the infer_ratio_d0_d1_unitary() routine among other utilities.

Classes¶

`KonCorrectionMLP`	G_genes MLPs indépendants entraînés en parallèle via joblib.
`GeneRegulatoryODE_softmax`	ODE model for gene regulatory dynamics with generalized softmax-based kon.

Functions¶

`build_kon_fn`(ks, theta_inter, bias[, device])	Return a function kon(X_numpy_or_torch) -> numpy array (batch, G)
`fit_scale_theta`(X_prot, kon_beta, bias, theta_inter, ...)	Find a single scale factor minimising the total MSE across all genes jointly.
`infer_ratio_d0_d1_unitary`(→ tuple[numpy.ndarray, ...)	Estimate ε_i = d1_i/d0_i from ODE residuals (bursty-PDMP variance matching).
`inference_degradation_prot`(→ tuple[numpy.ndarray, ...)	Estimate degradation rates and scaling factors from protein time-course data.
`infer_ratio_d0_d1_full`(→ tuple[numpy.ndarray, ...)	Infer d1/d0 from full simulated trajectories via regularised LS.
`predict_trajectory`(ode_func, X0, t_span[, method, ...])	Simulate a trajectory given an initial state and trained ODE model.
`compare_trajectories_umap`(ode_func, X_prot, times[, ...])	Compare real and simulated trajectories using UMAP projection.

Module Contents¶

CardamomOT.inference.degradations.build_kon_fn(ks, theta_inter, bias, device='cpu')¶: Return a function kon(X_numpy_or_torch) -> numpy array (batch, G) The function accepts either numpy arrays or torch tensors; it returns numpy.

class CardamomOT.inference.degradations.KonCorrectionMLP(G_tot: int, G_genes: int, hidden_dim: int = 8)¶

Bases: torch.nn.Module

G_genes MLPs indépendants entraînés en parallèle via joblib. Gène i prend [P_i, kon_beta_i] → delta_mu_i. At init: delta_mu = 1 everywhere.

class CardamomOT.inference.degradations.GeneRegulatoryODE_softmax(G, d_init, ks, theta_inter, bias, n_stimuli=1, stim_vals=None, device='cpu', kon_mlp=None, lambda_scale=1000.0)¶

Bases: torch.nn.Module

ODE model for gene regulatory dynamics with generalized softmax-based kon. Learns gene-specific degradation rates (d) and scale factors. Optionally applies a KonCorrectionMLP multiplicative correction (Harissa branch).

set_ratio_interpolation(t_start, t_end, ratio_t0, ratio_t1, lambda_mlp)¶

Configure time-interpolated ratio mix for the next odeint call.

Parameters:

t_start – boundaries of the current interval (floats).
t_end – boundaries of the current interval (floats).
ratio_t0 – mean observed g = kon_harissa/kon_beta at t_start / t_end, shape (G_genes,).
ratio_t1 – mean observed g = kon_harissa/kon_beta at t_start / t_end, shape (G_genes,).
lambda_mlp – mix weight (1 = pure interpolation, 0 = pure MLP).

forward(t, X)¶: Compute dX/dt for a given state X at time t. Includes learned scaling of theta_inter and bias. When kon_mlp is set, applies a frozen per-gene multiplicative correction. If set_ratio_interpolation has been called, the correction is mixed with the linearly interpolated observed training ratios at weight lambda_mlp.

CardamomOT.inference.degradations.fit_scale_theta(X_prot, kon_beta, bias, theta_inter, ks, ns, samples_data=None)¶

Find a single scale factor minimising the total MSE across all genes jointly.

Parameters:

X_prot – (N, G) protein levels.
kon_beta – (N, G) mixture-inferred kon values.
bias – (G, n_modes-1) or (n_samples, G, n_modes-1) GRN basal.
theta_inter – (G, G, n_modes-1) GRN interaction tensor.
ks – (n_modes, G) softmax burst-rate amplitudes.
ns – number of stimulus columns.
samples_data – (N,) per-cell sample index when bias is 3-D, else None.

Returns:

float, the jointly optimal scale (same for all genes).

Return type:

scale_theta

CardamomOT.inference.degradations.infer_ratio_d0_d1_unitary(X_prot, times, bias, theta_inter, ks, d_learned_temporal, k1_vec, n_stimuli=1, stim_schedule=None, samples_data=None, method='dopri5', rtol=1e-06, atol=1e-08, lambda_deg=0.0, prior_eps=None, outlier_quantile=0.95, min_kon=1e-06, eps_min=0.01, eps_max=100.0, verbose=True, kon_mlp=None, scale=1.0) → tuple[numpy.ndarray, numpy.ndarray]¶

Estimate ε_i = d1_i/d0_i from ODE residuals (bursty-PDMP variance matching).

Theoretical background.

In the bursty PDMP actually simulated (see BurstyPDMP), the burst parameters at interval cnt are:

burst rate : λ_i = k1_i · kon_norm_i(P) · d0_i
burst size : Exp(k1_i · d0_i / (scale · d1_i)) i.e. mean burst = scale · d1_i / (k1_i · d0_i)

This ensures the mean-field ODE limit is dP_i/dt = d1_i(kon_i/k1_i − P_i).

By Campbell’s theorem (leading term for small dt):

\[\operatorname{Var}[\Delta X_i] \approx \lambda_i \cdot E[B_i^2] \cdot dt = 2\,dt\cdot k_{on,i}(P)\cdot\mathrm{scale}\, \frac{d_{1,i}^2}{k_{1,i}^2\,d_{0,i}} = \varepsilon_i \cdot h_{ic}\]

where \(\varepsilon_i = d_{1,i}/d_{0,i}\) and

\[h_{ic} = 2\,dt\cdot\mathrm{scale}\cdot \frac{d_{1,i}}{k_{1,i}^2}\cdot k_{on,i}(P_c)\]

The regularised method-of-moments estimator is:

\[\varepsilon_i = \frac{\sum_c (X_{\mathrm{pred},c} - X_{1,c})^2 + \lambda\,\varepsilon_{\mathrm{prior},i}} {\sum_c h_{ic} + \lambda}\]

Parameters:

X_prot – (N, G) protein observations (all timepoints).
times – (N,) time label for each row.
bias – GRN bias — (T-1, G, n_nets) or (T-1, n_samples, G, n_nets).
theta_inter – (T-1, G, G, n_nets) GRN interaction tensor.
ks – (n_modes, G) softmax burst-rate amplitudes (already multiplied by scale_proteins).
d_learned_temporal – (T-1, G) learned protein degradation rates d1.
k1_vec – (G,) max burst rate × scale (= k1 * scale).
n_stimuli – Number of stimulus columns (default 1).
stim_schedule – {float: np.ndarray(ns)} stimulus values per time.
samples_data – (N,) per-cell sample index (optional).
method – ODE solver settings.
rtol – ODE solver settings.
atol – ODE solver settings.
lambda_deg – Tikhonov weight; 0 = pure MoM, large → prior.
prior_eps – (G,) prior for ε = d1/d0. Defaults to all-ones.
outlier_quantile – Fraction of cells to keep (sorted by residual).
min_kon – Minimum kon value for a cell to contribute.
eps_min – Output clipping bounds for ε = d1/d0.
eps_max – Output clipping bounds for ε = d1/d0.
verbose – Log per-interval diagnostics.
kon_mlp – Optional KonCorrectionMLP (Harissa branch).
scale – Protein scale (self.scale_proteins). Default 1.0.

Returns:

(T-1, G) per-interval ε = d1/d0 estimates. eps_global : (G,) global ε pooled over all intervals.

Return type:

eps_temporal

CardamomOT.inference.degradations.inference_degradation_prot(X_prot, times, bias, theta_inter, ks, d=None, n_epochs=500, lr=0.01, method='dopri5', rtol=1e-06, atol=1e-08, print_every=50, batch_size=None, verbose=True, n_stimuli=1, stim_schedule=None, scale_proteins=1.0, samples_data=None, kon_mlp=None, lambda_scale=1000.0, lambda_deg=0.0, lambda_mlp=0.0, g_obs_train=None) → tuple[numpy.ndarray, numpy.ndarray]¶

Estimate degradation rates and scaling factors from protein time-course data.

When bias is 3-D (n_samples, G, n_modes-1) and samples_data is provided, one ODE module is created per sample with its own bias while d_param and scale_param are shared, so a single optimizer refines the shared kinetics from all samples jointly.

Returns:: Learned degradation rates, shape (G,). scale_learned: Learned scaling parameters, shape (G,).
Return type:: d_learned

CardamomOT.inference.degradations.infer_ratio_d0_d1_full(X_prot, times, bias, theta_inter, ks, d_learned, k1_vec, kon_mlp, prior_d1d0=None, n_stimuli=1, stim_schedule=None, samples_data=None, method='dopri5', rtol=1e-05, atol=1e-07, n_steps=10, lambda_deg=0.0, lambda_mlp=0.5, g_obs_train=None, min_h=0.0001, verbose=True, clip_lo=0.01, clip_hi=100.0) → tuple[numpy.ndarray, numpy.ndarray]¶

Infer d1/d0 from full simulated trajectories via regularised LS.

Theoretical background.

To first order in d1/d0, the MLP correction factor satisfies:

\[1 - g_i(P(t)) \approx \frac{d_{1,i}}{d_{0,i}} \cdot h_i(P(t))\]

where g_i = kon_harissa_i / kon_beta_i (learned by kon_mlp) and

\[h_i(P) = \sum_j \frac{\partial\ln kon_i}{\partial X_j} \cdot (kon^{\rm eff}_j - X_j)\]

where the d_{1,j} factor is not included in h_i (uniform-d1 approximation: d_{1,j} \approx d_{1,i} for all j, so it becomes part of the ratio d_{1,i}/d_{0,i} that the LS estimates directly).

Fitting strategy.

Instead of evaluating only at the terminal state, this function simulates full trajectories with n_steps intermediate evaluation points and fits r_i = d1_i/d0_i by minimising over all (cell, timestep) pairs:

\[\mathcal{L}(r_i) = \sum_{c,t} \bigl[y_{ict} - r_i\,h_{ict}\bigr]^2 + \lambda\,(r_i - r_{\rm prior,i})^2\]

with y_{ict} = 1 - g_i(P_c(t)). The closed-form solution is:

\[r_i^* = \frac{\sum h\,y + \lambda\,r_{\rm prior}}{\sum h^2 + \lambda}\]

Only (cell, timestep) pairs where h_i and y_i have the same sign (i.e. where the formula predicts r_i > 0) and where |h_i| > min_h are included.

Usage in base.py:

ratios_temporal, ratios_global = infer_ratio_d0_d1_full(...)
# temporal (one ratio per interval):
self.ratios[cnt] = 1.0 / ratios_temporal[cnt]
# global (one shared ratio):
self.ratios[:] = (1.0 / ratios_global)[None, :]

Parameters:

X_prot – (N, G) protein observations across all timepoints.
times – (N,) time label for each observation.
bias – GRN bias — (G, n_modes-1) or (n_samples, G, n_modes-1).
theta_inter – GRN interactions — (G, G, n_modes-1) or (n_samples, G, G, n_modes-1).
ks – Mode amplitudes, shape (n_modes, G).
d_learned – Protein degradation rates, shape (G,) (constant) or (T-1, G) (temporal). When 2-D, interval cnt uses d_learned[cnt] for the ODE simulation.
k1_vec – Kept for interface compatibility (not used internally).
kon_mlp – Trained KonCorrectionMLP.
prior_d1d0 – (G,) prior for d1/d0 used when lambda_deg > 0. Defaults to all-ones (quasi-stationary limit).
n_stimuli – Number of stimulus columns (default 1).
stim_schedule – Dict {float_time: np.ndarray(ns,)}; defaults to all-ones when None.
samples_data – (N,) per-cell sample index.
method – ODE solver method (default "dopri5").
rtol – ODE solver tolerances.
atol – ODE solver tolerances.
n_steps – Number of intermediate ODE evaluation points per interval (default 10; total = n_steps+1 incl. t0, t1).
lambda_deg – Tikhonov regularisation weight — same role as in inference_degradation_prot(). 0 = pure LS; large values → prior.
lambda_mlp –
Mix weight between the linearly interpolated observed training ratios and the MLP prediction. For a simulated state at time t in [t1, t2] the effective ratio used as y_i = 1 - g_eff is:

\[g_{\rm eff} = \lambda_{\rm mlp} \Bigl[ \frac{t_2-t}{t_2-t_1}\,r_1 + \frac{t-t_1}{t_2-t_1}\,r_2 \Bigr] + (1-\lambda_{\rm mlp})\,g(P,\,k_{\rm on}(P))\]

where r_1 / r_2 are the mean observed MLP ratios at the training timepoints t_1 / t_2 (shape (G_genes,)). 1.0 = pure linear interpolation of observed g; 0.0 = pure MLP prediction (original behaviour).
min_h – Minimum |h_i| to include a (cell, timestep) pair.
verbose – Log per-interval diagnostics.
clip_lo – Final clamp on the estimated d1/d0.
clip_hi – Final clamp on the estimated d1/d0.

Returns:

(T-1, G) per-interval d1/d0 estimates. ratios_global : (G,) global d1/d0 pooled over all intervals.

Stimulus columns ([:ns]) are fixed to 1.0 in both outputs.

Return type:

ratios_temporal

CardamomOT.inference.degradations.predict_trajectory(ode_func, X0, t_span, method='dopri5', rtol=1e-06, atol=1e-08, stim_vals=None)¶: Simulate a trajectory given an initial state and trained ODE model.

CardamomOT.inference.degradations.compare_trajectories_umap(ode_func, X_prot, times, method='dopri5')¶: Compare real and simulated trajectories using UMAP projection.