wnetdeconv

Spectral deconvolution via Wasserstein optimal transport.

Given an empirical spectrum and a library of theoretical component spectra, wnetdeconv finds the mixture proportions that minimise the total Wasserstein transport cost between the empirical signal and the weighted sum of components. The inner problem at each set of proportions is solved exactly as a min-cost flow (via pylmcf / LEMON), giving an exact piecewise-linear objective with exact gradients — suitable for gradient- based outer optimisation with scipy.

Supports 1-D spectra (NMR chemical shift, m/z) and higher-dimensional data (e.g. m/z + retention time).

Installation

pip install wnetdeconv

Dependencies: pylmcf, wnet, numpy, scipy. Optional: pyopenms for loading featureXML files.

Concepts

Spectra as distributions

A spectrum is a set of (position, intensity) pairs. In 1-D (NMR chemical shift, m/z) use Spectrum_1D; for higher-dimensional data (m/z + retention time) use Spectrum with a (d, n) positions array.

from wnetdeconv import Spectrum_1D

empirical = Spectrum_1D([1.0, 2.0, 3.0], [10.0, 25.0, 15.0])
component = Spectrum_1D([1.0, 2.0, 3.0], [1.0, 2.0, 1.0])

Transport cost

Matching a unit of intensity from an empirical peak at position p to a theoretical peak at position q costs distance(p, q). Peaks that cannot be matched cheaply are instead routed to a trash node at a fixed penalty.

max_distance caps the farthest match considered; anything farther is cheaper to trash. trash_cost (or the asymmetric pair experimental_trash_cost / theoretical_trash_cost) sets that penalty.

Scaling

Internally all intensities and costs are scaled to integers for the MCF solver. Intensities are quantised automatically by a p95-quantile policy (wnet.scaling.WNetDeconvScaler); cost quantisation is chosen by the network itself against its integer budget. Construction raises ValueError if quantisation would silently discard more than 20% of any spectrum's total intensity (see allow_intensity_loss). Pass scale_factor to override the automatic scaling explicitly. The scipy ftol stop criterion is derived from the actual scale factors, so the outer loop stops as soon as further improvement is below the resolution the integer network can deliver.

The former precision parameter is deprecated and has no effect; passing a non-default value raises a DeprecationWarning.

Solvers

DeconvSolver — unconstrained baseline

Solves the network at a given point and exposes total_cost() and gradient(). Optimisation (via optimize(), L-BFGS-B) minimises cost with only non-negativity bounds.

from wnetdeconv import DeconvSolver, Spectrum_1D
from wnet.distances import DistanceMetric

emp  = Spectrum_1D([1.0, 100.0], [10.0, 30.0])
t1   = Spectrum_1D([1.0],        [2.0])   # optimal proportion: 5
t2   = Spectrum_1D([100.0],      [3.0])   # optimal proportion: 10

solver = DeconvSolver(
    empirical_spectrum=emp,
    theoretical_spectra=[t1, t2],
    distance=DistanceMetric.LINF,
    max_distance=10.0,
    trash_cost=100.0,
)

result = solver.optimize()
print(result.x)   # ~[5. 10.] (exact up to intensity quantization)

You can also drive the solver manually — useful when embedding it in your own optimisation loop:

solver.set_point([5.0, 10.0])
print(solver.total_cost())   # 0.0
print(solver.gradient())     # [200. 300.]

Note that the gradient is not zero at the optimum: the objective is piecewise linear, and at a kink gradient() reports the right-derivative (here, the cost of trashing the surplus theoretical mass that one more unit of each proportion would add).

ConstrainedSolver — total-mass equality

Adds the constraint Σ wₛ · Iₛ = I_emp so that the mixture exactly accounts for all empirical intensity. Uses SLSQP. Drop-in replacement for DeconvSolver; call optimize() the same way.

from wnetdeconv import ConstrainedSolver

solver = ConstrainedSolver(
    empirical_spectrum=emp,
    theoretical_spectra=[t1, t2],
    distance=DistanceMetric.LINF,
    max_distance=10.0,
    trash_cost=100.0,
)
result = solver.optimize()

optimize_cutting_plane() — kink-safe alternative optimizer

The objective is convex piecewise linear in the proportions, and descent methods (L-BFGS-B, SLSQP) can stall on the kinks between linear pieces when several components' spectra overlap heavily. optimize_cutting_plane() (on DeconvSolver and ConstrainedSolver) runs Kelley's cutting-plane method instead: each evaluation contributes a supporting plane and the next iterate minimizes the accumulated plane model over the feasible polytope (a small LP per iteration). It returns the best evaluated point together with a model lower bound and gap; the mass-balance constraint is carried natively by the LP. Typically converges in 10-40 evaluations and never returns a point with higher cost than it evaluated. By default (polish=True) the answer is then refined by the class's descent optimizer warm-started from the cutting-plane point, and the cheaper of the two is returned — the plane model navigates the kinks, the descent finishes the final linear piece, and the comparison makes the combination monotone (never worse than either stage alone).

Key parameters

Distance metrics

From wnet.distances.DistanceMetric:

• L1 — sum of absolute coordinate differences (Manhattan / taxicab)
• L2 — Euclidean distance
• LINF — maximum absolute coordinate difference (Chebyshev); dual of the W₁ earth-mover distance used by masserstein

Loading MS data (featureXML)

from wnetdeconv import Spectrum

emp = Spectrum.FromFeatureXML("sample.featureXML")   # requires pyopenms

Architecture

wnetdeconv
├── Spectrum / Spectrum_1D   — data containers (extend wnet.Distribution)
├── DeconvSolver             — core: builds WassersteinNetwork, exposes cost + gradient
├── ConstrainedSolver        — adds total-mass equality, uses SLSQP
│   └── MagnetsteinSolver    — magnetstein-style: all spectra normalised to sum 1;
│                              with MTD_th, defaults to independent two-sided trash
│                              (dualdeconv3/4 semantics; independent_trash=False
│                              restores the annihilating model)
├── MassersteinSolver2       — reproduces masserstein's dualdeconv2 (one-sided trash)
└── MassersteinSolver4       — reproduces dualdeconv4 (independent two-sided trash);
                               requires wnet >= 1.3.0 (independent trash)

The underlying min-cost flow is provided by wnet (network construction) and pylmcf (LEMON-based MCF algorithms, including warm-restart Network Simplex).

Related work

Sister paper: WNetAlign applies the same truncated Wasserstein / network simplex machinery to MS and NMR spectral alignment:

Król J, Bochenek M, Jopa S, Kazimierczuk K, Gambin A, Startek MP (2026). WNetAlign: fast and accurate spectra alignment using truncated Wasserstein distance and network simplex. Briefings in Bioinformatics, 27(3), bbag247. https://doi.org/10.1093/bib/bbag247

License

MIT