Python api update

README.md

@@ -9,7 +9,17 @@ for iterative fitting against experimental data.
 
 ### Installation
 
-Requires a Rust toolchain and Python ≥ 3.9.
+Install straight from GitHub with [uv](https://docs.astral.sh/uv/). The
+extension is built from source, so a
+[Rust toolchain](https://www.rust-lang.org/tools/install) and Python ≥ 3.9 are
+required:
+
+```sh
+uv pip install "git+https://github.com/mlund/pripps.git#subdirectory=pripps-py"
+```
+
+To also get the example structures and form-factor assets used below, clone the
+repo and install in editable mode instead:
 
 ```sh
 git clone https://github.com/mlund/pripps.git
@@ -29,79 +39,90 @@ model = pr.multipole_model(
     qmin=0.0,
     qmax=0.5,
     nq=500,
-    excluded="foxs",         			   # FoXS-style excluded volume (16π Fraser)
-    volume_scale=1.0,        			   # fit knob c₁ (scales excluded volume)
-    hydration="sasa",  # FoXS-style SASA-weighted hydration shell
+    excluded="foxs",          # or pripps.Excluded.Foxs — 16π Fraser excluded volume
+    hydration="sasa",         # or pripps.Hydration.Sasa
     probe=1.8,  # SASA probe radius (Å)
-    contrast_density=0.03,  # fit knob c₂ (scales hydration)
     bulk_electron_density=0.334,  # e/ų, pure water
 )
 
-# Returns a dict of numpy arrays
-out = model.intensity(volume_scale=1.0, contrast_density=0.03)
-print(out["q"][:3], out["total"][:3])
+# Fit knobs (volume_scale, contrast_density) are applied here, not at build time.
+res = model.intensity(volume_scale=1.0, contrast_density=0.03)   # -> pripps.Intensity
+print(res.q[:3], res.total[:3])
 ```
 
-The returned dict contains:
-
-| Key | Description |
-|-----|-------------|
-| `q` | Scattering vector magnitudes (Å⁻¹) |
-| `total` | Total intensity I(q) |
-| `atoms` | Protein-only (vacuum) contribution |
-| `excluded` | Excluded-volume contribution |
-| `hydration` | Hydration-shell contribution |
+`intensity()` returns a typed `Intensity` with NumPy-array attributes:
 
-### Fitting against experimental data
+| Attribute | Description |
+|-----------|-------------|
+| `.q` | Scattering vector magnitudes (Å⁻¹) |
+| `.total` | Total intensity I(q) |
+| `.atoms` | Protein-only (vacuum) contribution |
+| `.excluded` | Excluded-volume contribution |
+| `.hydration` | Hydration-shell contribution |
 
-`multipole_model` does the expensive per-structure setup once and
-caches the intermediate scattering amplitudes. Subsequent
-`model.intensity(volume_scale, contrast_density)` calls only apply
-the fit knobs to the cached amplitudes — typically under a
-millisecond each — so wrapping them in `scipy.minimize` or a grid
-search is fluent.
-
-For small structures or exact reference checks, `debye_model` accepts
-the same solvent keyword arguments and exposes the same
-`model.intensity(volume_scale, contrast_density)` interface. It caches
-the exact Debye self and cross terms for atoms, excluded volume, and
-hydration before fitting.
-`direct_model` does the same for the explicit discrete-q transform,
-caching the complex amplitudes before recombining and averaging
-duplicate q magnitudes.
-
-A full working example is in `pripps-py/examples/fit_lysozyme.py`,
-which fits lysozyme against the shipped `tests/lysozyme/lyzexp.dat`
-experimental dataset and reproduces FoXS's χ²/N ≈ 0.20:
+Cached models always populate `.atoms` / `.excluded` / `.hydration` (zero
+where a contribution is disabled). They are `None` only when there is no
+per-contribution breakdown — currently trajectory averages (the one-shot
+`Pripps.intensity(..., xtc=...)`).
 
-```sh
-python pripps-py/examples/fit_lysozyme.py
-```
+### Fitting against experimental data
 
-The core pattern:
+`multipole_model` does the expensive per-structure setup once and caches the
+intermediate amplitudes; `model.intensity(volume_scale, contrast_density)` then
+only applies the fit knobs — typically under a millisecond. There are two ways
+to fit:
 
 ```python
 import numpy as np
-from scipy.optimize import minimize
 
 # experimental data: q (Å⁻¹), I(q), σ(q)
-data = np.loadtxt("experimental.dat", comments="#")
-q_exp, i_exp, sigma = data[:, 0], data[:, 1], data[:, 2]
+q_exp, i_exp, sigma = np.loadtxt("experimental.dat", comments="#", unpack=True)
+
+# (a) built-in fit — returns a typed FitResult:
+fit = model.fit(q_exp, i_exp, sigma)
+print(fit.volume_scale, fit.contrast_density, fit.chi2, fit.scale)
 
+# (b) manual fit with scipy for full control:
+from scipy.optimize import minimize
 def chi2(params):
-    c1, c2 = params
-    out = model.intensity(volume_scale=c1, contrast_density=c2)
-    theory = np.interp(q_exp, out["q"], out["total"])
+    out = model.intensity(*params)
+    theory = np.interp(q_exp, out.q, out.total)
     scale = np.sum(theory * i_exp / sigma**2) / np.sum((theory / sigma) ** 2)
     return np.sum(((scale * theory - i_exp) / sigma) ** 2) / len(q_exp)
+result = minimize(chi2, x0=[1.0, 0.03],
+                  bounds=[(0.95, 1.12), (-2.0, 4.0)], method="L-BFGS-B")
+```
+
+`debye_model` (exact, for small structures) and `direct_model` (explicit
+discrete-q) accept the same solvent arguments and expose the same
+`intensity` / `fit` / `q_values` interface.
 
-result = minimize(
-    chi2,
-    x0=[1.0, 0.2],
-    bounds=[(0.95, 1.12), (-2.0, 4.0)],
-    method="L-BFGS-B",
+A full working example is in `pripps-py/examples/fit_lysozyme.py`, which fits
+lysozyme against the shipped `tests/lysozyme/lyzexp.dat` and reproduces FoXS's
+χ²/N ≈ 0.20 (demonstrating both fitting styles):
+
+```sh
+python pripps-py/examples/fit_lysozyme.py
+```
+
+### Trajectory averaging
+
+To average I(q) over an MD trajectory, use the one-shot `Pripps.intensity`
+(the cached models are single-structure). The `model=` argument selects the
+calculator; trajectory averages carry no per-contribution breakdown:
+
+```python
+avg = pr.intensity(
+    structure_file="protein.pdb",   # topology / template frame
+    model="multipole",              # or "debye" / "direct", or a pripps.Scheme
+    xtc="trajectory.xtc",
+    stride=10,                      # use every 10th frame
+    weights="weights.dat",          # optional per-frame weights
+    excluded="foxs",
+    hydration="sasa",
+    contrast_density=0.03,
 )
-print(f"c1={result.x[0]:.4f}, c2={result.x[1]:.4f}, χ²/N={result.fun:.4f}")
+print(avg.q, avg.total)
 ```
 
 ### Solvent models
@@ -110,9 +131,9 @@ print(f"c1={result.x[0]:.4f}, c2={result.x[1]:.4f}, χ²/N={result.fun:.4f}")
 |-----------|--------|-------------|
 | `excluded` | `"disabled"`, `"fraser"`, `"foxs"`, `"grid"`, `"voronoi"`, `"per-kind"` | Excluded-volume model |
 | `hydration` | `"disabled"`, `"sasa"`, `"grid"`, `"voronoi"` | Hydration-shell model |
-| `volume_scale` | float (default 1.0) | Scales excluded volume (PepsiSAXS `r₀`) |
+| `volume_scale` | float (default 1.0) | Scales excluded volume (PepsiSAXS `r₀`). In Python this is a fit knob passed to `model.intensity()` / `model.fit()`, not the model builder |
 | `excluded_spacing` | float (default 1.0) | Grid spacing (Å) for `excluded="grid"` |
-| `contrast_density` | float (default 0.03) | Fractional hydration-layer excess density (c₂ in `A − c₁·C + c₂·B`) — pepsi's `d_rho / ρ_water`, ~0.1 for a 10 % excess layer |
+| `contrast_density` | float (default 0.03) | Fractional hydration-layer excess density (c₂ in `A − c₁·C + c₂·B`) — pepsi's `d_rho / ρ_water`, ~0.1 for a 10 % excess layer. In Python a fit knob on `model.intensity()` / `model.fit()` |
 | `probe` | float (default: 1.8 for SASA/grid, 6.0 for voronoi) | Water probe radius in Å |
 | `spacing` | float (default 3.10516) | Grid spacing (Å) for `hydration="grid"`; the bulk-water linear spacing, so each cell holds ~1 water molecule |
 | `thickness` | float (default 3.0) | Shell thickness (Å) for `hydration="grid"` |

assets/stovgaard2010.py

@@ -10,6 +10,6 @@
     if residue == "q" or residue == "Unnamed: 0":
         continue
     print(f"{residue}: {{ formfactor: !Table [", end="")
-    for q, I in zip(df["q"], df[residue]):
-        print(f"[{q}, {I}], ", end="")
+    for q, intensity in zip(df["q"], df[residue]):
+        print(f"[{q}, {intensity}], ", end="")
     print("]}\n")

pripps-py/examples/batch_single_bead_fit_sasbdb.py

@@ -13,8 +13,8 @@
 as testing_single-bead_approximation.ipynb:
   - atomfile: assets/poly_amino_acid.yaml (or poly_martini.yaml for Martini3)
   - excluded: per-kind using the atomfile's excluded_solvent entries
-  - hydration: disabled
-  - optimizer: L-BFGS-B with bounds [(0.95, 1.12), (-2.0, 4.0)] or change them in the script
+  - hydration: grid
+  - fit: built-in `model.fit(q, I, sigma)` (pripps' library optimizer)
 """
 
 from __future__ import annotations
@@ -25,7 +25,6 @@
 
 import numpy as np
 import pripps
-from scipy.optimize import minimize
 
 
 def load_exp_data(path: Path, qmax: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
@@ -54,12 +53,6 @@ def load_exp_data(path: Path, qmax: float) -> tuple[np.ndarray, np.ndarray, np.n
     return data[:, 0], data[:, 1], data[:, 2]
 
 
-def weighted_scale(theory: np.ndarray, i_exp: np.ndarray, sigma: np.ndarray) -> float:
-    num = np.sum(theory * i_exp / sigma**2)
-    den = np.sum((theory / sigma) ** 2)
-    return float(num / den)
-
-
 def fit_single(
     pr: pripps.Pripps,
     xyz_path: Path,
@@ -76,38 +69,25 @@ def fit_single(
         nq=nq,
         excluded="per-kind",
         hydration="grid",
-        volume_scale=1.0,
         probe=1.8,
-        contrast_density=0.03,
         bulk_electron_density=0.334,
-        spacing = 3.8,
+        spacing=3.8,
     )
 
-    def chi2(params: np.ndarray) -> float:
-        c1, c2 = float(params[0]), float(params[1])
-        out = model.intensity(volume_scale=c1, contrast_density=c2)
-        theory = np.interp(q_exp, out["q"], out["total"])
-        scale = weighted_scale(theory, i_exp, sigma)
-        return float(np.sum(((scale * theory - i_exp) / sigma) ** 2) / len(q_exp))
-
-    result = minimize(
-        chi2,
-        x0=[1.0, 0.2],
-        bounds=[(-2.0, 3.0), (-2.0, 4.0)],
-        method="L-BFGS-B",
-    )
+    # Built-in fit of (volume_scale, contrast_density) against the data.
+    fit = model.fit(q_exp, i_exp, sigma)
 
-    c1_opt, c2_opt = float(result.x[0]), float(result.x[1])
-    out = model.intensity(volume_scale=c1_opt, contrast_density=c2_opt)
-    theory_interp = np.interp(q_exp, out["q"], out["total"])
-    scale = weighted_scale(theory_interp, i_exp, sigma)
-    i_fit = scale * theory_interp
+    out = model.intensity(
+        volume_scale=fit.volume_scale, contrast_density=fit.contrast_density
+    )
+    theory_interp = np.interp(q_exp, out.q, out.total)
+    i_fit = fit.scale * theory_interp
 
     return {
-        "c1": c1_opt,
-        "c2": c2_opt,
-        "chi2": float(result.fun),
-        "scale": scale,
+        "c1": fit.volume_scale,
+        "c2": fit.contrast_density,
+        "chi2": fit.chi2,
+        "scale": fit.scale,
         "q_exp": q_exp,
         "i_exp": i_exp,
         "sigma": sigma,

pripps-py/examples/fit_lysozyme.py

@@ -5,6 +5,11 @@
 excluded-volume formula (16π denominator), and SASA-weighted
 hydration. Reproduces FoXS's χ²/N ≈ 0.20 on lysozyme.
 
+Demonstrates BOTH fitting styles:
+  (a) the built-in `model.fit(...)` convenience, and
+  (b) a manual `scipy.optimize.minimize` over the typed `Intensity`
+      result (`.q` / `.total`) for full control.
+
 Run from the repository root (after `cd pripps-py && uv pip install -e .`):
 
     python pripps-py/examples/fit_lysozyme.py
@@ -32,60 +37,67 @@
 
 # Prepare the Multipole model — expensive work runs once here.
 # Uses FoXS-compatible Cromer-Mann form factors with the 16π
-# Fraser excluded-volume formula and SASA hydration.
+# Fraser excluded-volume formula and SASA hydration. Model choice
+# may be a string ("foxs") or an enum (pripps.Excluded.Foxs).
 pr = pripps.Pripps(atomfile="assets/foxs_formfactors.yaml")
 model = pr.multipole_model(
     structure_file="tests/lysozyme/6lyz.pdb",
     qmin=0.0,
     qmax=float(q_exp[-1]),
     nq=500,
-    excluded="foxs",
-    hydration="sasa",
+    excluded=pripps.Excluded.Foxs,
+    hydration=pripps.Hydration.Sasa,
 )
 
+# (a) Built-in fit — one call, returns a typed FitResult.
+fit = model.fit(q_exp, i_exp, sigma)
+print("Built-in model.fit():")
+print(f"  {fit!r}")
+print("  FoXS reference: χ²/N ≈ 0.20")
+
 
+# (b) Manual fit with scipy for full control over the objective.
+# The typed Intensity gives numpy `.q` / `.total` directly.
 def chi2(params):
-    """χ² between the theoretical curve and experimental data."""
     c1, c2 = params
     out = model.intensity(volume_scale=c1, contrast_density=c2)
-    theory = np.interp(q_exp, out["q"], out["total"])
+    theory = np.interp(q_exp, out.q, out.total)
     scale = np.sum(theory * i_exp / sigma**2) / np.sum((theory / sigma) ** 2)
     return np.sum(((scale * theory - i_exp) / sigma) ** 2) / len(q_exp)
 
 
-# Grid search to find a good starting point
-print("Grid search over (volume_scale, contrast_density)...")
-best_chi2 = float("inf")
-best_params = (1.0, 0.0)
-for c1 in np.linspace(0.95, 1.12, 18):
-    for c2 in np.linspace(-2.0, 4.0, 61):
-        val = chi2([c1, c2])
-        if val < best_chi2:
-            best_chi2 = val
-            best_params = (c1, c2)
-print(f"  Best grid point: c1={best_params[0]:.3f}, c2={best_params[1]:.3f}, χ²/N={best_chi2:.4f}")
-
-# Refine with L-BFGS-B
 result = minimize(
     chi2,
-    x0=list(best_params),
+    x0=[1.0, 0.0],  # neutral start — independent of (a), to confirm they agree
     bounds=[(0.95, 1.12), (-2.0, 4.0)],
     method="L-BFGS-B",
 )
-c1_opt, c2_opt = result.x
-print(f"  Optimised:       c1={c1_opt:.4f}, c2={c2_opt:.4f}, χ²/N={result.fun:.4f}")
-print(f"  FoXS reference:  χ²/N ≈ 0.20")
-
-# Write the fitted curve
-out = model.intensity(volume_scale=c1_opt, contrast_density=c2_opt)
-theory = np.interp(q_exp, out["q"], out["total"])
-scale = np.sum(theory * i_exp / sigma**2) / np.sum((theory / sigma) ** 2)
+print("Manual scipy.minimize():")
+print(f"  c1={result.x[0]:.4f}, c2={result.x[1]:.4f}, χ²/N={result.fun:.4f}")
 
+# Write the fitted curve (using the built-in fit's parameters).
+out = model.intensity(
+    volume_scale=fit.volume_scale, contrast_density=fit.contrast_density
+)
+theory = np.interp(q_exp, out.q, out.total)
 np.savetxt(
     "fit_lysozyme.csv",
-    np.column_stack([q_exp, i_exp, scale * theory]),
+    np.column_stack([q_exp, i_exp, fit.scale * theory]),
     delimiter=",",
     header="q,i_exp,i_fit",
     comments="",
 )
-print(f"  Fitted curve written to fit_lysozyme.csv (scale={scale:.6e})")
+print(f"Fitted curve written to fit_lysozyme.csv (scale={fit.scale:.6e})")
+
+# --- Trajectory averaging (one-shot; the only path that reads an XTC) ---
+# Uncomment with a real trajectory to average I(q) over frames:
+# avg = pr.intensity(
+#     structure_file="tests/lysozyme/6lyz.pdb",
+#     model="multipole",
+#     xtc="trajectory.xtc",
+#     stride=10,
+#     excluded="foxs",
+#     hydration="sasa",
+#     contrast_density=0.03,
+# )
+# print(avg.q, avg.total)   # note: trajectory averages have no per-contribution breakdown