Struct JointPoissonObjective 

pub struct JointPoissonObjective<'a> {
    pub model: &'a dyn FitModel,
    pub o: &'a [f64],
    pub s: &'a [f64],
    pub c: f64,
    pub active_mask: Option<&'a [bool]>,
    pub open_background: Option<&'a [f64]>,
    pub sample_background: Option<&'a [f64]>,
}

Joint-Poisson objective.

Wraps an effective count-ratio FitModel (which produces T_eff,i = model.evaluate(θ)) together with the observed open-beam counts O_i, sample counts S_i, and proton-charge ratio c = Q_s / Q_ob.

Physical count-response contract

This is a low-level API: it cannot inspect a dyn FitModel to determine how instrument resolution was applied. The supplied model must already return the physically valid effective sample/open count response on the observed bins. With an instrument response operator R and incident spectrum Φ, that response is

T_eff = R[Φ · T] / R[Φ].

A post-hoc broadened transmission R[T] is not this response and must never be supplied as though it were. In particular, do not directly wrap a resolution-bearing TransmissionFitModel, because that model returns R[T]. A no-resolution transmission model remains valid because R is then the identity; a custom model that already returns the exact effective ratio is also valid. Otherwise, callers must model the open and sample response arms separately or use the guarded pipeline entry point, which rejects unsupported counts-plus-resolution combinations.

The caller is responsible for ensuring o, s, and model.evaluate() output all have the same length.

Fields

model: &'a dyn FitModel

Effective count-ratio model: evaluate(θ) → T_eff(E).

See the struct-level physical count-response contract. In particular, this must not be a post-hoc broadened transmission R[T].

o: &'a [f64]

Open-beam counts per bin.

s: &'a [f64]

Sample counts per bin.

c: f64

Proton-charge ratio c = Q_s / Q_ob. Must be strictly positive.

active_mask: Option<&'a [bool]>

Optional per-bin active mask (SAMMY EMIN/EMAX-equivalent fit-energy-range restriction). When Some(m), only bins where m[i] is true contribute to the deviance / gradient / Fisher information; the model is still evaluated on the full grid so resolution broadening at the boundaries is correct. When None, all bins are active (default behaviour).

Length must equal o.len(); the GUI / pipeline dispatch builds it from the configured [E_min, E_max] against the energy grid.

open_background: Option<&'a [f64]>

Expected background counts in the OPEN-beam acquisition, per bin.

None is the background-free case and keeps the binomial reduction documented above. Some switches to the Poisson form, because with a background in either arm the flux no longer cancels out of the conditional likelihood.

These are counts for the acquisition they belong to, already binned — a measured dark or blocked-beam reference after its own run normalization, not a transmission-level curve.

sample_background: Option<&'a [f64]>

Expected background counts in the SAMPLE acquisition, per bin.

Distinct from Self::open_background: the sample scatters neutrons and emits gammas, so equal backgrounds is not a physical case. Their difference is exactly what the two arms can separate.

Implementations

impl<'a> JointPoissonObjective<'a>

pub fn n_data(&self) -> usize

Number of data bins.

pub fn n_active(&self) -> usize

Number of active data bins — n_data when no mask is set, or the count of true entries in active_mask otherwise.

pub fn n_informative(&self) -> usize

Number of informative active bins: active bins with a nonzero count total O_i + S_i > 0. A zero-total bin is degenerate under the conditional-binomial model — its profiled rate is zero and it contributes exactly zero deviance for every parameter value — so counting it as a degree of freedom deflates deviance_per_dof (and the opt-in scale_by_chi2 σ inflation) by the empty-bin fraction. The exact detector-time route makes wide acquisition windows with many empty bins routine, so deviance-per-dof reporting must use THIS count.

pub fn profile_lambda(&self, t_i: f64, o_i: f64, s_i: f64) -> f64

Closed-form profile MLE for the per-bin flux: λ̂ = c·(O+S) / (1+c·T).

Guards: when 1 + c·T ≤ ε, returns 0 to avoid division blow-up.

pub fn profile_lambda_with_background( &self, t_i: f64, o_i: f64, s_i: f64, b_o: f64, b_s: f64, ) -> f64

Profile MLE for the per-bin flux WITH backgrounds present.

With E[O] = λ/c + b_o and E[S] = λ·T + b_s, the score equation

-1/c + O/(λ + c·b_o) - T + S·T/(λ·T + b_s) = 0

clears to a quadratic a λ² + b λ + d = 0 with

a = T·(1/c + T)
b = (1/c + T)·(b_s + c·b_o·T) - T·(O + S)
d = c·b_o·b_s·(1/c + T) - O·b_s - S·T·c·b_o

a > 0 and d <= 0 whenever the counts are non-negative, so there is exactly one non-negative root and it is the + branch. Setting b_o = b_s = 0 gives d = 0 and λ̂ = -b/a = c(O+S)/(1+cT), the background-free closed form — which is the check that this reduces correctly, and is asserted by zero_backgrounds_reproduce_the_binomial_profile.

Solved with the numerically stable quadratic form: the + branch of (-b + sqrt(disc))/(2a) cancels catastrophically when b > 0, so it is evaluated as -2d/(b + sqrt(disc)) there.

pub fn profile_lambda_per_bin( &self, t: &[f64], ) -> Result<Vec<f64>, FittingError>

Vector form of profile_lambda.

Validates t.len() == o.len() == s.len() and c > 0; returns FittingError::LengthMismatch / InvalidConfig rather than the previous .zip() truncate-and-pretend behaviour (which would silently shrink the output to min(t.len(), o.len(), s.len())).

pub fn deviance_from_transmission(&self, t: &[f64]) -> Result<f64, FittingError>

Conditional binomial deviance at the given transmission vector.

D = 2 · Σ [ S·ln(S/(Np)) + O·ln(O/(N(1−p))) ] with p = cT/(1+cT), N = O+S, and x·ln(x/0) → 0.

Near invalid or numerically tiny transmission values, the per-bin evaluation (binomial_deviance_term) uses t.max(POISSON_EPSILON) to clamp T away from zero before entering the logarithms and the 1/(1+cT) factor. This avoids singular logs and division-by-zero but is a piecewise clamp, not a smooth quadratic extrapolation — D(T) is C⁰ at the clamp boundary, not C¹. In practice this is adequate because the optimizer’s transmission values come from a FitModel that keeps T bounded well above POISSON_EPSILON for physically plausible density / nuisance parameter values.

pub fn deviance(&self, params: &[f64]) -> Result<f64, FittingError>

Evaluate the deviance at parameter vector θ by calling the model.

pub fn deviance_gradient_analytical( &self, params: &[f64], free_param_indices: &[usize], ) -> Result<Option<Vec<f64>>, FittingError>

Analytical gradient of the deviance w.r.t. the free parameters.

Returns None if the transmission model does not provide an analytical Jacobian — callers should fall back to deviance_gradient_fd.

Gradient derivation: with p_i = cT_i/(1+cT_i) and N_i = O_i+S_i,

d D / d T_i = −2 · (S_i − O_i·c·T_i) / (T_i · (1 + c·T_i))

then chain-rule with the transmission Jacobian J_{i,j} = ∂T_i / ∂θ_{f(j)} where f(j) is the j-th free parameter index.

pub fn fisher_information( &self, params: &[f64], free_param_indices: &[usize], ) -> Result<Option<FlatMatrix>, FittingError>

Fisher information for free parameters (Gauss-Newton curvature of D).

Uses the expected-info form

h_i ≡ ∂² D / ∂ T_i² ≈ 2 · (O_i + S_i) · c / (T_i · (1 + c·T_i)²)

(derived from logit-link binomial Var(S|N) = N p (1−p) and d logit(p) / dT = 1/T, scaled by 2 since D = −2 L). Then

I(θ){j,k} = Σ_i h_i · J{i,j} · J_{i,k}.

Returns None if the transmission model does not provide an analytical Jacobian.

pub fn fisher_information_fd( &self, params: &mut ParameterSet, fd_step: f64, ) -> Result<Option<FlatMatrix>, FittingError>

Finite-difference Fisher information.

Fallback for callers whose transmission model does not implement FitModel::analytical_jacobian — i.e., when Self::fisher_information would return None. Builds the transmission Jacobian column-by-column via central differences and assembles

I(θ)_{j,k} = Σ_i h_i · J_{i,j} · J_{i,k}

where h_i = ∂² D / ∂ T_i² is the per-bin deviance curvature 2·(O_i + S_i)·c / (T_i·(1 + c·T_i)²) (Fisher-scoring form derived from binomial logit-link Var(S | N) = N·p·(1−p) with d logit p / dT = 1/T — see the module-level docstring §Model). Returns Ok(None) only if the base model evaluation itself fails.

pub fn deviance_gradient_fd( &self, params: &mut ParameterSet, fd_step: f64, ) -> Result<Vec<f64>, FittingError>

Finite-difference gradient of the deviance.

Central differences on each free parameter. Used as a fallback when the model has no analytical Jacobian. params is a mutable ParameterSet so we can respect bounds via clamp().