Struct ExpGaussian 

pub struct ExpGaussian;

Exponentially-modified Gaussian (EMG) — a Gaussian convolved with a one-sided exponential, the canonical asymmetric/tailing chromatography & spectroscopy peak.

A · (γ/2) · exp[γ(c−x) + (γσ)²/2] · erfc[(c + γσ² − x)/(σ√2)]

Parameters (in order): [amplitude, center, sigma, gamma]

• gamma is the exponential decay rate of the tail (toward high x); as γ→0 the shape approaches a Gaussian.

Numerical stability (no clamp)

The naive form computes exp(arg_exp)·erfc(z), which overflows to inf·0 → NaN for arg_exp > 709 (e.g. γσ > 37). Instead we use the algebraic identity arg_exp − z² = −(x−c)²/(2σ²) and split on the sign of z:

• z ≥ 0: A·(γ/2)·exp(−(x−c)²/(2σ²))·erfcx(z) — both factors are bounded (erfcx(z) ∈ (0,1], the Gaussian ≤ 1), so there is no overflow.
• z < 0: A·(γ/2)·exp(arg_exp)·erfc(z) — here arg_exp < 0, so the exp cannot overflow, and erfc(z) ∈ (1,2).

The two branches are continuous at z = 0 (erfcx(0) = erfc(0) = 1 and the Gaussian factor equals exp(arg_exp) there). A final is_finite guard remains as belt-and-suspenders. The numpy benchmark oracle uses the identical split with scipy.special.erfcx, so numpy↔Rust parity holds to machine precision.

Trait Implementations

impl Model for ExpGaussian

fn jacobian(&self, x: &[f64], params: &[f64]) -> Vec<f64>

Analytical Jacobian of the EMG (exponentially-modified Gaussian).

Define: e_arg = γ(c−x) + (γσ)²/2, u = (c + γσ² − x)/(σ√2), E = exp(e_arg), C = erfc(u), g_u = (2/√π)·exp(−u²).

Then f = A·(γ/2)·E·C.

∂f/∂A = (γ/2)·E·C ∂f/∂center = A·(γ/2)·E·[ γ·C − g_u/(σ√2) ] ∂f/∂sigma = A·(γ/2)·E·[ γ²σ·C − g_u·(γ√2 − u/σ) ] ∂f/∂gamma = (A/2)·E·C + A·(γ/2)·E·[ (c−x+γσ²)·C − g_u·σ/√2 ]

When the overflow-clamped region returns f=0, all derivatives are 0.

Numerical stability

The Jacobian uses the same identity as eval: E·C = gauss·erfcx(u) for z ≥ 0 (overflow-free) E·C = exp(e_arg)·erfc(u) for z < 0 (safe, e_arg < 0 here)

In both branches, E·g_u = gauss·(2/√π) where gauss = exp(−(x−c)²/(2σ²)), because the identity e_arg − u² = −(x−c)²/(2σ²) holds exactly.

fn eval(&self, x: &[f64], params: &[f64]) -> f64

Evaluate the model at coordinate x with the given params.

fn jacobian_into(&self, x: &[f64], params: &[f64], out: &mut [f64])

Fill a pre-allocated slice with Jacobian values (one entry per parameter).

fn param_names(&self) -> Vec<Cow<'static, str>>

Ordered parameter names, matching the layout expected by eval and jacobian.

fn eval_slice_into(&self, xs: &[f64], params: &[f64], out: &mut [f64])

Batch evaluation of a 1-D model: fill out[i] = eval([xs[i]], params).

fn jac_slice_into(&self, xs: &[f64], params: &[f64], out: &mut [f64])

Batch Jacobian for a 1-D model: fill out in row-major layout [i * params.len() + j] = d(model)/d(params[j]) at xs[i].

fn n_dims(&self) -> usize

Number of coordinate dimensions consumed from x. Defaults to 1.