Feature request: composable floor-normalization + analytic per-series scaling for relative / arbitrary-unit data

[wshlavacek]

Summary

Two small, composable per-series normalization primitives that would let PyBNF express the
log / relative objectives common to arbitrary-unit systems-biology data using standard tokens,
instead of a bespoke objective class:

1. Floor normalization — x' = x + rho * max(x) per series (rho a small constant, e.g. 0.03).
2. Analytic per-series scaling — profile out the optimal per-series multiplicative scale
   (a.k.a. hierarchical / optimal scaling), rather than fitting one scale parameter per series.

Together with the existing lognormal objective, these compose the "sum of squared log-differences
of geometric-mean-normalized trajectories" objective used across the Lipniacki-lab NF-kB models —
and, in various forms, most arbitrary-unit fluorescence/blot fits.

Motivation (worked example)

Jaruszewicz-Blonska et al., PLoS ONE 2023 18(6):e0286416 fit a reduced NF-kB model to reproduce
the original Lipniacki-2004 model by minimizing (their Eq 7):

J = sum_i ( ln y_i^R  -  ln y_i^O )^2

where each measured series is normalized (Methods) by

x'_i = x_i + 0.03 * max(x)     # a measurement-noise floor
y_i  = x'_i / geomean(x')       # arbitrary units -> divide by the geometric mean

The geomean division makes the objective invariant to an overall multiplicative scale per series
(the reduced model is non-dimensional, the original is dimensional); the floor keeps the log finite
where a series is legitimately zero (e.g. active IKK at rest).

We reconstructed this as two edition-2 jobs (BNGL-Models pybnf-jobs/Jaruszewicz-Blonska-2023/).
The closest we can get with current tokens is objective = norm_sos + normalization = peak + a
rho-floored target — a first-order analog that reproduces the paper's reported AMD exactly (Gate 3a)
and recovers 12/13 parameters (Gate 3b), but leaves the paper's least-identifiable parameter
(epsilon) a few-fold looser than it should be under the exact objective, because the analog's
minimum is not at the same point.

Current surface and the gaps

• Log residuals already exist: objective = lognormal = Gaussian on log10 (objective.py:1186,
  1461). (Minor: it desugars to sigma = read_exp_file _SD, so a fixed-sigma log objective
  currently needs _SD columns or a noise_model line; a lognormal shortcut with sigma = fix_at 1
  would help this compose.)
• Normalization families are init | peak | zero | unit (data.py:127, 401-514) — none adds a
  max-fraction floor, and none profiles out a per-series scale.
• The only user-defined-objective hook is analytical / model-free. objective = expression | callable (ADR-0050; objective.py:1502-1514 -> DirectPassObjective; config.py:978+ synthesizes
  an ExpressionModel/CallableModel with no .bngl and no simulation) optimizes a closed form over
  the parameters — it cannot see a simulated trajectory, so it cannot express Eq 7. There is therefore
  no way to express a custom sim-vs-data objective short of a new @register_objfunc class. Hence
  this request is for the two composable primitives rather than a paper-specific objective.

Proposed features

1. Floor normalization

A per-series transform x' = x + rho * max(x) (default rho = 0.03): a well-defined log/relative
objective on data that legitimately touches zero, plus down-weighting of near-zero measurement noise.
Suggested surface, composable with the existing normalizations:

normalization <obs> = floor 0.03
normalization <obs> = floor 0.03, peak      # floor, then peak

Design note: this must be applied identically to sim and data. Today normalization is sim-only
(the .exp is pre-normalized by the user); a floor is only meaningful applied symmetrically to both.

2. Analytic per-series scaling (optimal / hierarchical scaling)

For each series (experiment x observable), analytically profile out the optimal multiplicative scale
before scoring, so an overall scale difference between model and data is not penalized (arbitrary
units). Closed forms: linear least squares c* = sum(s*d)/sum(s*s); log-space
c* = geomean(d)/geomean(s) (= mean-centering the logs — the paper's geomean case). This is the
"hierarchical optimization" / analytic-scaling approach (Weber et al. 2011; Loos et al. 2018;
implemented in pyPESTO / AMICI / dMod), and it avoids one free scale parameter per series (which
would be 54 extra parameters in the combination job). Suggested surface:

normalization <obs> = scale                 # analytic optimal scale (family-appropriate)

or an observable-level observable: X, scale: analytic.

Composition

With both, the paper's exact objective is spelled with standard tokens:

objective = lognormal                        # fixed sigma
normalization <obs> = floor 0.03, scale

applied to both sim and data — no bespoke objective class, and reusable for any arbitrary-unit
relative-data fit.

Notes

• Both are per-series (per experiment x observable) transforms. The scale case maps naturally to
  PEtab estimated observableParameters / hierarchical scaling; the floor is a preprocessing step —
  worth stating each one's PEtab-export status when implemented.
• Could be split into two issues; filed together because they share the relative-data-log-objective
  motivation and compose for it.