Metadata-Version: 2.4
Name: wnetdeconv
Version: 0.11.0
Summary: Python implementation of spectral deconvolution using Wasserstein metric
Author-Email: Michał Startek <michal.startek@mimuw.edu.pl>
Maintainer-Email: Michał Startek <michal.startek@mimuw.edu.pl>
License-Expression: MIT
License-File: LICENCE
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Mathematics
Classifier: Topic :: Software Development :: Libraries :: Python Modules
Classifier: Development Status :: 4 - Beta
Project-URL: Homepage, https://github.com/michalsta/wnetdeconv
Project-URL: Repository, https://github.com/michalsta/wnetdeconv.git
Requires-Python: >=3.10
Requires-Dist: pylmcf>=1.2.1
Requires-Dist: numpy
Requires-Dist: scipy
Requires-Dist: wnet>=1.3.0
Requires-Dist: nanobind-backend>=1.0
Provides-Extra: extras
Requires-Dist: pyopenms; extra == "extras"
Provides-Extra: pytest
Requires-Dist: pytest; extra == "pytest"
Requires-Dist: pandas; extra == "pytest"
Description-Content-Type: text/markdown
Dynamic: Requires-Dist

# 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](https://github.com/michalsta/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

```bash
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.

```python
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.

```python
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:

```python
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.

```python
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

| Parameter | Applies to | Description |
|---|---|---|
| `max_distance` | all | Maximum peak-to-peak match distance. Also sets the sparsity of the internal network in 1-D. |
| `trash_cost` | all | Symmetric penalty for unmatched peaks. |
| `experimental_trash_cost` | `DeconvSolver` | Per-unit penalty for discarding empirical mass. |
| `theoretical_trash_cost` | `DeconvSolver` | Per-unit penalty for discarding theoretical mass. |
| `precision` | all | Deprecated, no effect; non-default values raise `DeprecationWarning`. |
| `scale_factor` | all | Override automatic intensity/cost quantisation. |
| `smooth_sigma` | `MagnetsteinSolver` | Gaussian pre-smoothing width in position units (default 0 = off). On evenly gridded profile 1H NMR, ~0.02 ppm (a typical linewidth) consistently helps when mixture and reference lineshapes differ (e.g. field-inhomogeneity distortion); much wider blurs small components away. Leave off for centroided data. |

## 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)

```python
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](https://github.com/michalsta/wnet) (network construction) and
[pylmcf](https://github.com/michalsta/pylmcf) (LEMON-based MCF algorithms,
including warm-restart Network Simplex).

## Related work

**Sister paper:** [WNetAlign](https://github.com/michalsta/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
