95 add bortfeld model and wilkens models

[witNie]

This pull request introduces several improvements to the pyamtrack project, focusing on expanding functionality, clarifying Python API conventions, and providing practical usage examples. The most significant changes are the addition of a new module to the build process, detailed documentation of Python argument naming conventions, and a new example notebook demonstrating the dose_bortfeld function.

Build system update:

• Added proton_models to the PYAMTRACK_TARGETS in CMakeLists.txt, ensuring that the proton models module is built and linked as part of the project.

Documentation enhancements:

• Introduced a comprehensive naming_convention.md document detailing the rules and current state of Python argument names and calling conventions for the nanobind-based API. This includes general naming rules, argument types, default values, validation, and a dictionary of all public arguments.

Examples and usage:

• Added a new Jupyter notebook examples/dose_bortfeld.ipynb that demonstrates how to use the proton_models.dose_bortfeld function, including parameter setup, function call, and plotting the resulting dose profile.

[grzanka]

🤖 AI-generated explainer (model: Claude Opus 4.8), prepared to give the team shared vocabulary and physics context for this PR.
Sources: T. Bortfeld, An analytical approximation of the Bragg curve for therapeutic proton beams, Med. Phys. 24, 2024 (1997) — the single-Bragg-peak model that dose_bortfeld actually implements; and J. J. Wilkens & U. Oelfke, Analytical linear energy transfer calculations for proton therapy, Med. Phys. 30, 806 (2003) — the LET model. Please sanity-check against the physics before treating any of it as gospel.

---

Proton dose & LET models — a primer for the team

0. Why these models exist (60-second version)

We treat tumours by firing protons into the body. The one fact that makes protons special:

A proton deposits most of its energy right before it stops, at a well-defined depth, and almost nothing beyond. That burst is the Bragg peak.

X-rays dump most of their energy near the surface and keep going; protons let you place the "explosion" inside the tumour and spare tissue behind it.

The two models here answer two questions:

1. Bortfeld → How much energy lands at each depth? (the dose curve)
2. Wilkens → How concentrated is that energy at each depth? (the LET curve)

Both are analytical: closed-form formulas evaluated directly, instead of a slow Monte-Carlo simulation. They run in microseconds — which is what makes them usable inside an optimisation loop.

---

1. Core vocabulary (read once, refer back forever)

Term	Plain meaning	Unit	In code
Depth	how far into the material, along the beam	cm	z_cm, depth_cm
Fluence	how many particles crossed a 1 cm² window	cm⁻²	fluence_cm2
Dose	energy absorbed per kg of tissue — "how hard you got hit"	Gray, Gy (J/kg)	return of dose_bortfeld
Stopping power S	energy a single proton loses per unit distance, −dE/dx	MeV/cm	internal
Range R	distance a proton of a given energy travels before stopping	cm	internal
Residual range r	distance a proton still has left	cm	internal (Wilkens)
Bragg peak	the dose spike near the end of the range	—	the shape dose_bortfeld makes
Range straggling	protons don't all stop at exactly the same depth; the peak blurs	cm (a σ)	fed by sigma_E_MeV
LET	Linear Energy Transfer — local energy density along the track	keV/µm	return of the Wilkens fn
RBE	Relative Biological Effectiveness — biological damage per Gy	—	(motivation only)

Two unit traps:

• Dose is Gy, LET is keV/µm. Different things, different units — don't mix them.
• keV/µm = (MeV/cm) × 0.1. Wilkens computes in MeV/cm internally, then ×0.1.

---

2. The single assumption everything is built on: the range–energy power law

Every formula in both papers falls out of one empirical rule:

R = α · E^p

In words: a proton's range is its starting energy raised to a power.

• E = starting kinetic energy (MeV), R = how deep it gets (cm)
• α, p = two fitted constants. For protons in water, p ≈ 1.77, α ≈ 0.0022.

Worked example — a 150 MeV proton in water:

R = 0.0022 × 150^1.77 ≈ 0.0022 × 7100 ≈ 15.6 cm

Measured value ≈ 15.8 cm. A two-constant formula lands within a couple of millimetres — that's why the whole approach works (a log–log plot of range vs energy is a straight line).

For programmers: p is a modelling knob, not a law of nature. Its "correct" physical value is ~1.8, but a fit may use a different value to quietly compensate for ignored effects. If you see a magic 1.77/1.8 in the code, that's this exponent.

---

3. The Bortfeld single Bragg curve — what dose_bortfeld computes

Take the range–energy law and turn the crank:

1. A proton at depth d has just enough energy left to cover the remaining distance R − d. Invert the power law → energy at every depth.
2. Dose ∝ stopping power = −dE/dx. Differentiate:

D(d) ∝ 1 / (R − d)^(1 − 1/p)

Look at that denominator: as depth d → range R, (R − d) → 0, so dose shoots up. That spike is the Bragg peak, from one line of algebra.

Mathematically it's a singularity (infinite dose) — unphysical. Two real effects round it into the measured peak, and these are exactly the extra parameters in the code:

• Range straggling (sigma_E_MeV) — protons stop at slightly different depths; modelled as a Gaussian blur. This is why the function needs an energy spread, not just an energy.
• Nuclear interactions / tail (eps) — a small fraction of protons hit nuclei early and never reach the peak. This adds a little dose to the entrance region and shaves the peak. eps is that fraction (default ≈ 0.03, ~3%).

Reading the signature with new eyes

dose = pm.dose_bortfeld(
    z_cm,         # depth(s) where we want the dose          [cm]
    fluence_cm2,  # how many protons per cm²                 [1/cm²]
    E_MeV,        # beam energy → sets WHERE the peak lands   [MeV]
    sigma_E_MeV,  # energy spread → how BLURRY the peak is    [MeV]
    material,     # 1 = liquid water                          [id]
    eps,          # nuclear tail fraction (~0.03)             [—]
)

Mental model: energy moves the peak, spread softens it, eps adds the entrance tail, fluence scales the whole thing up and down.

On the eps = -1.0 review comment

-1 is a sentinel value: "I didn't specify — use your default." libamtrack sees the negative number and substitutes 0.03. It's a C idiom (no Optional/None in plain C). The reviewer's point is that on the Python side we can expose the real default (0.03) directly, so users see a meaningful number instead of a magic -1. Both behave identically; one is self-documenting. Same story for sigma_E_MeV = -1 → 0.01·E.

---

4. Beam geometry: pencil beam vs. large field — does it matter?

Yes, and it's worth being explicit, because it defines what fluence_cm2 and the returned dose actually mean.

Both models are central-axis, 1-D depth models: dose (Bortfeld) and LET (Wilkens) as a function of depth only. Wilkens says it outright — "central axis of broad proton beams." Neither describes the lateral (off-axis) profile; that would be a separate Gaussian you convolve in.

They assume lateral charged-particle equilibrium: at the scoring point, as many protons scatter out sideways as scatter in. That holds in two practically equivalent situations:

1. Large homogeneous field (e.g. passively scattered), measured with a small detector on the central axis, away from the field edge. Lateral scatter cancels, so on-axis dose depends only on depth.
2. A single pencil beam measured by a chamber wider than the beam — i.e. the integrated / laterally-summed depth dose (IDD), as from a large Bragg-peak chamber. Integrating over the lateral plane recovers the same curve, because the protons that scatter off-axis are still caught by the wide chamber.

These two give the same depth-dose shape — the broad-beam ≈ integrated-pencil equivalence (superposition, assuming the lateral profile just shifts with depth and lateral nuclear transport is small).

Where it does NOT apply: the on-axis depth dose of a narrow pencil beam read by a small detector. There, lateral scattering spreads the beam with depth and bleeds dose off the axis, so the measured central-axis curve falls faster than these models predict. That case needs the full 2-D/3-D pencil-beam model (depth model × lateral Gaussian).

For our wrappers: treat fluence_cm2 as the planar fluence of a broad parallel beam (equivalently, the per-pencil fluence integrated laterally), and the output as central-axis-of-large-field / integrated depth dose. It is not the on-axis dose of a thin pencil read by a small detector. Same caveat for Wilkens LET.

---

5. LET — the concept behind the review comment

Dose tells you how much energy landed. LET tells you how concentrated it was along each track. Biology cares about concentration, not just total.

Analogy: two people each drink 10 units of alcohol (same "dose") — one over a week, one in a night (different "LET"). Same total, very different effect. Densely-deposited energy (high LET) causes clustered, hard-to-repair DNA damage.

Formally, LET = stopping power, in keV/µm:

• Entrance: fast protons, energy spread thinly → low LET (~0.5 keV/µm).
• Near the peak: protons crawling, dumping energy in a tiny space → high LET (several keV/µm).

That's why the tail end of a proton beam is biologically more potent than dose alone suggests — and why we compute LET at all (it feeds RBE).

The crux: there is no single "LET" at a point — there are two averages

At any depth you don't have one proton, you have a spectrum of them (different energies from straggling and the beam's initial energy width), each with its own stopping power. So "the LET here" needs an average over the spectrum, and there are two standard ones:

Quantity	Averaged by	Symbol
Track-averaged LET	particle count (fluence) — plain mean of S	LET_t (Lt)
Dose-averaged LET	dose contribution — slow, high-S protons count more	LET_d (Ld)

LET_t = ∫ w(r)·S(r) dr  / ∫ w(r) dr           (mean of S)
LET_d = ∫ w(r)·S²(r) dr / ∫ w(r)·S(r) dr       (S-weighted mean of S)

w(r) is the spectrum (how many protons have residual range r). LET_d has S² on top, so it up-weights the slow high-stopping-power protons. Therefore:

LET_d ≥ LET_t always, and the gap grows near the Bragg peak. For monoenergetic protons the two are equal.

In radiobiology, dose-averaged LET (LET_d) is the one people actually use for RBE.

⛔ Why the review comment is correct

The PR currently names the function let_dose and documents it as "the LET dose." That phrase doesn't exist in the field and reads like "LET × dose," which is a different quantity. The C function it calls is AT_LET_d_Wilkens_keV_um_* — the _d means dose-averaged LET.

The fix is purely naming/wording:

Current (wrong)	Correct
let_dose	let_dose_averaged (or dose_averaged_let)
"the LET dose"	"the dose-averaged LET (LET_d), in keV/µm"

If the track-averaged variant gets wrapped later, it'd be let_track_averaged → AT_LET_t_Wilkens_*. Keeping _d / _t visible in the Python name mirrors the C API and removes the ambiguity for good.

---

6. The Wilkens model — LET, built directly on Bortfeld

Wilkens' insight (2003): if Bortfeld gets the dose curve from the power-law range–energy rule, the same machinery gives the LET curve. He reuses:

• the same power law R = α E^p (p = 1.77),
• the same stopping-power fit (S(r) ∝ r^(1/p − 1)),
• a Gaussian proton spectrum at each depth (mean residual range R₀ − z, width σ from straggling + initial beam energy spread).

Plug those into the two averaging integrals and they solve in closed form (using parabolic cylinder functions — the LET analogue of Bortfeld's arctan/log for dose). Feed in beam energy + energy spread, get LET_t and LET_d at any depth, validated against Monte-Carlo to within 0.5 keV/µm.

One detail you'll see referenced as R (the regularization, 2 mm): the raw stopping-power formula blows up as r → 0 (a nearly-stopped proton has "infinite" stopping power on paper). Unphysical, so Wilkens averages S over the last 2 mm of range to tame it — the LET equivalent of how sigma/eps tame Bortfeld's dose singularity.

let_d = pm.let_dose_averaged(   # suggested name
    depth_cm,          # depth(s)                  [cm]
    material,          # 1 = water                 [id]
    energy_MeV,        # beam energy               [MeV]
    energy_spread_MeV, # initial energy width σ_E  [MeV]
)                      # → dose-averaged LET       [keV/µm]

---

7. How the pieces fit together

        range–energy power law   R = α·E^p          ← the single shared assumption
                     │
        ┌────────────┴─────────────┐
        ▼                           ▼
  Bortfeld dose                Wilkens LET
  (energy WHERE?)              (energy how CONCENTRATED?)
        │                           │
  single Bragg peak           LET_t  (fluence-weighted mean of S)
  (Bortfeld 1997 →            LET_d  (dose-weighted mean of S)   ← what we wrap
   dose_bortfeld)               ↑ "dose-averaged LET", NOT "LET dose"

  both: central-axis / broad-beam (large field on axis,
        or integrated pencil-beam depth dose) — not narrow-pencil-on-a-small-detector

---

8. Concrete numbers to sanity-check the code

For a 150 MeV proton beam in water (material = 1):

Quantity	Expect roughly	Why
Peak depth (range)	~15.6 cm	0.0022 × 150^1.77
Entrance dose	small, flat	far from the peak
Peak dose	several × entrance	the 1/(R−d) blow-up, rounded by sigma/eps
LET_d in entrance	~0.5 keV/µm	fast protons, dilute
LET_d at distal edge	several keV/µm, rising sharply	protons crawling to a stop
LET_d vs LET_t	LET_d ≥ LET_t	the S² weighting

If dose_bortfeld puts the peak near 15–16 cm and the LET curve rises steeply at the distal edge while staying ~0.5 keV/µm in the entrance, the wrappers are behaving.

Minor aside spotted in let_example.ipynb: pm.let_dose(depths, 1, 150, 1.5) raises std::bad_cast when depths is built via list(numpy_array) — the C++ side handles a Python list, a float, or an nb::ndarray<double>, but a list of numpy.float64 matches none cleanly. Separate from the vocabulary fixes, but worth a follow-up.

---

Generated by Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[Copilot]

Pull request overview

Copilot reviewed 13 out of 13 changed files in this pull request and generated 1 comment.

Suppressed comments (7)

examples/let_wilkens.ipynb:144

• cartesian=True is not a valid keyword argument. The bound parameter is cartesian_product. Passing cartesian=True raises TypeError for an unexpected keyword. (This line also uses the non-existent particle_models module and the wrong positional argument order noted above.)

        "results = pyamtrack.particle_models.let_wilkens(depths, pyamtrack.materials.water_liquid, 150, 0.01, averaging=pyamtrack.particle_models.Averaging.TRACK, cartesian=True)"

examples/dose_bortfeld.ipynb:27

• This example notebook was committed with a stale failed execution: the first cell's saved output is a ModuleNotFoundError: No module named 'matplotlib', and all subsequent cells have no outputs. As shipped, this example does not demonstrate the Bortfeld model. Please re-run it in an environment with matplotlib installed (or clear the outputs) before committing.

        {
          "output_type": "error",
          "ename": "ModuleNotFoundError",
          "evalue": "No module named 'matplotlib'",
          "traceback": [
            "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
            "\u001b[31mModuleNotFoundError\u001b[39m                       Traceback (most recent call last)",
            "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[1]\u001b[39m\u001b[32m, line 2\u001b[39m\n\u001b[32m      1\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m numpy \u001b[38;5;28;01mas\u001b[39;00m np\n\u001b[32m----> \u001b[39m\u001b[32m2\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m matplotlib.pyplot \u001b[38;5;28;01mas\u001b[39;00m plt\n\u001b[32m      3\u001b[39m \n\u001b[32m      4\u001b[39m \u001b[38;5;28;01mimport\u001b[39;00m pyamtrack\n\u001b[32m      5\u001b[39m \n",
            "\u001b[31mModuleNotFoundError\u001b[39m: No module named 'matplotlib'"
          ]
        }
      ],

src/proton_models/proton_models_bindings.cpp:116

• Typo in comment: "constucted" should be "constructed".

	// enum constucted by chaining

src/proton_models/proton_models_bindings.cpp:58

• The Raises section refers to the internal C++ names z_cm, E_MeV, and sigma_E_fraction, but the Python-facing parameters (documented above) are depth_cm, energy_MeV, and energy_spread_fraction. Use the public names here so the user-facing docstring is self-consistent.

            If z_cm < 0, E_MeV outside [0.1, 10000.0], sigma_E_fraction outside (0, 1),
            eps outside [0, 1), or a material ID is not a known material.

src/proton_models/proton_models_bindings.cpp:14

• The new proton_models module (dose_bortfeld here and let_wilkens below) has no automated tests, while the other bound modules are covered (e.g. tests/test_stopping.py, tests/test_materials.py). Consider adding a tests/test_proton_models.py covering scalar/list/array inputs, the cartesian_product path, the averaging selector, and the validation error cases (negative depth, out-of-range energy/spread/eps, invalid material).

  m.def("dose_bortfeld", &dose_bortfeld, nb::arg("depth_cm"), nb::arg("fluence_cm2"), nb::arg("energy_MeV"),

src/proton_models/dose_bortfeld.h:22

• The hardcoded range [1, 24] is misleading. Material validation is table-driven (validate_material_argument), and materials_bindings.cpp now enumerates actual IDs via get_ids() rather than assuming contiguous 1..N, so material IDs are not guaranteed to be [1, 24]. Reference the canonical source instead.

 * - material: int material_no in [1, 24] OR pyamtrack.materials.Material object

CMakeLists.txt:154

• These two top-level commands were indented by 2 spaces, which is inconsistent with the surrounding top-level statements in this file (e.g. lines 156-157). Remove the indentation to keep the file consistent.

  # Create the Python module from the source file
  nanobind_add_module(_core src/main.cpp)

[grzanka]

Critical review of #203 — Bortfeld dose & Wilkens LET models

How this was checked: I built the head commit (6ae2cc9) into a wheel on Linux/CPython 3.11 (GSL 2.x, libamtrack master), installed it, ran the existing suite (87 passed, 1 skipped — no regressions), and then probed the new API and cross-checked the physics against libamtrack/library/src/AT_ProtonAnalyticalModels.c and AT_DataMaterial.{c,h}. Every claim below is either a code citation or a reproduced observation. Numbers quoted are from that build.

What works well

• Wheel builds and imports cleanly; adding proton_models to PYAMTRACK_TARGETS is all that was needed (the foreach/file(GLOB ...) wiring is generic) ✅
• The physics looks right. 150 MeV in water: Bragg peak at z = 15.4 cm, just proximal to the model's range α·E^p = 15.69 cm (α = 0.00231, p = 1.761 from the libamtrack table); entrance LET_d ≈ 0.54 keV/µm; LET_d = 3.17 vs LET_t = 2.30 keV/µm at 15 cm. All as expected.
• Hiding libamtrack's negative-sentinel convention (sigma_E_MeV < 0 → 0.01·E, eps < 0 → 0.03) behind real Python defaults is the right call and resolves @grzanka's earlier -1 comment.
• averaging="dose"/"track" correctly resolves @grzanka's "LET dose" comment — the naming now matches AT_LET_d/_t_Wilkens_*.
• materials_bindings.cpp switching from Material(i + 1) to actual get_ids() is a genuine robustness fix (it no longer assumes contiguous IDs).
• docs/naming_convention.md is a genuinely valuable document — clear, opinionated, and the right level of detail. My criticisms below are that this PR does not follow it.

Already flagged by the Copilot review and I agree, so I won't re-litigate: broken let_wilkens.ipynb calls, stale/failed notebook output, constucted typo, Raises section using internal names, missing tests, hardcoded [1, 24], CMakeLists.txt indentation. Details below where I have something to add.

---

Blocking

B1. The documented ValueError contract is wrong: list/array inputs raise RuntimeError

Both docstrings promise ValueError for out-of-range inputs. That holds only for scalars:

pm.let_wilkens(-1.0, 150.0)             # ValueError: depth_cm must be >= 0 ...          ✅
pm.let_wilkens([-1.0], 150.0)           # RuntimeError: Error processing 1-D NumPy array: ... ❌
pm.let_wilkens(np.array([-1.0]), 150.0) # RuntimeError ❌
pm.let_wilkens([-1.0], 150.0, cartesian_product=True)  # ValueError ✅
pm.dose_bortfeld([1.0], 1e8, 150.0, 0.01, 1, [5.0])    # RuntimeError (eps out of range) ❌

Cause: all validation lives inside the per-element lambda, and wrap_multiargument_function wraps the element loop in catch (const std::exception& e) { throw std::runtime_error(...); } (src/wrapper/multi_argument.h:126). std::invalid_argument — and nb::value_error, which derives from std::runtime_error — get swallowed and re-thrown as RuntimeError. The cartesian path has no such catch, hence the inconsistency between the two paths of the same function.

This is a pre-existing wrapper limitation, but #203 is the first code to rely on per-element validation, so it is this PR that makes the contract observable and wrong. Any test written from the docstring (pytest.raises(ValueError)) fails.

The over-broad catch also leaks results on every validation error (new double[input_length] at line 106; the owning capsule is only created at line 130). Both are fixed together:

// src/wrapper/multi_argument.h
std::unique_ptr<double[]> results(new double[input_length]);
try {
  for (size_t i = 0; i < input_length; i++) { /* ... */ results[i] = func(arguments_vector); }
} catch (const nb::cast_error&) {
  throw nb::type_error("1-D NumPy array dtype cannot be cast to double or input is not suitable.");
}
// No generic catch: let nanobind translate std::invalid_argument -> ValueError itself.
double* raw = results.release();
nb::capsule owner(raw, [](void* p) noexcept { delete[] (double*) p; });
return nb::ndarray<double, nb::numpy>(raw, {input_length}, owner).cast();

wrap_cartesian_product_function has the same leak (line 150) without the exception bug — worth the same unique_ptr treatment while you're in there.

B2. No tests for ~230 lines of new binding code

tests/ is untouched. The two new functions, the averaging selector, Averaging, valid_material_ids, the new material validation, and the changed Material constructor are all uncovered — and CI is green precisely because nothing exercises them. docs/naming_convention.md added in this PR says (step 6): "Add a test for the positional call and, when declared, the keyword call." Please hold the PR to its own rule.

A tests/test_proton_models.py covering the invariants I checked by hand would have caught B1, H1 and H4:

import numpy as np, pytest, pyamtrack
pm = pyamtrack.proton_models

def test_scalar_and_sequence_agree():
    assert pm.dose_bortfeld([10.0], 1e8, 150.0, 0.01)[0] == pytest.approx(
        pm.dose_bortfeld(10.0, 1e8, 150.0, 0.01))

def test_bragg_peak_near_range():
    z = np.linspace(0, 20, 2001)
    assert 15.0 < z[pm.dose_bortfeld(z, 1e8, 150.0, 0.01).argmax()] < 16.0

def test_let_d_at_least_let_t():
    z = np.linspace(0, 15, 151)
    lt = pm.let_wilkens(z, 150.0, averaging="track")
    ld = pm.let_wilkens(z, 150.0, averaging="dose")
    assert np.all(ld >= lt * (1 - 1e-8))          # see M6 for why the tolerance is needed

def test_material_id_and_object_agree():
    assert pm.let_wilkens(5.0, 150.0, 0.01, 1) == pm.let_wilkens(
        5.0, 150.0, 0.01, pyamtrack.materials.water_liquid)

@pytest.mark.parametrize("depth", [-1.0, [-1.0], np.array([-1.0])])
def test_negative_depth_is_value_error(depth):          # currently FAILS for list/array
    with pytest.raises(ValueError):
        pm.let_wilkens(depth, 150.0)

@pytest.mark.parametrize("material", [0, 9999, 1.0, "water"])
def test_invalid_material_rejected(material):
    with pytest.raises((ValueError, TypeError)):
        pm.let_wilkens(1.0, 150.0, 0.01, material)

def test_cartesian_shape():
    assert pm.let_wilkens([1.0, 2.0, 3.0], [150.0, 160.0],
                          cartesian_product=True).shape == (3, 2)

B3. examples/let_wilkens.ipynb cannot run, and isn't mentioned in the PR description

Confirmed against the built module — three independent breakages, in every let_wilkens cell:

pyamtrack.particle_models                       # AttributeError: module is proton_models
let_wilkens(depths, 1, 150, 0.01)               # TypeError: signature is (depth_cm, energy_MeV,
                                                #   energy_spread_fraction, material, ...)
let_wilkens(..., cartesian=True)                # TypeError: kwarg is cartesian_product

It reads as a leftover from an earlier API iteration. The PR description documents only the CMake change, the naming doc, and dose_bortfeld.ipynb — let_wilkens, the notebook, and the materials changes (which include a behaviour change, see H2) go unmentioned. Since the description becomes the squash-merge message, please bring it up to date. Also note the notebook's np.linspace(0, 20, 500) grid runs past the beam range, where let_wilkens returns NaN (H5).

---

High

H1. energy_spread_fraction must be strictly > 0 — a monoenergetic beam is impossible

pm.let_wilkens(1.0, 150.0, 0.0)     # ValueError: must be in range (0, 1)
pm.let_wilkens(1.0, 150.0, 1e-12)   # 0.5584433714241268  — accepted
pm.let_wilkens(1.0, 150.0, 0.01)    # 0.5584433714241268  — identical to 15 digits

Rejecting 0.0 while accepting 1e-12 is arbitrary, and zero spread is not a numerical hazard upstream. In AT_dose_Bortfeld_Gy_single the widths add in quadrature:

const double sigma_mono_cm = 0.012 * pow(range_cm, 0.935);  // range straggling, always > 0
const double sigma_cm = sqrt(sigma_mono_cm*sigma_mono_cm + pow(tmp_sigma_E_MeV*alpha*p*pow(E_MeV,p-1), 2));

so sigma_E_MeV = 0 yields sigma_cm = sigma_mono_cm > 0 — that is exactly Bortfeld's monoenergetic case, a physically meaningful and numerically safe limit that a user might reasonably want. Suggest relaxing to [0, 1) in both models and updating docstrings plus naming_convention.md. (Keep rejecting negatives — that is what stops users tripping the upstream sentinel.)

H2. Two different material-validity rules, one of them coupled to upstream table order

src/materials/materials.h states the intent: "This is the single point of material validation: callers must not re-check IDs themselves", then implements it as

if (id < 1 || AT_material_index_from_material_number(id) < 0)   // validate_material_argument

but src/materials/materials.cpp:8 was changed to a different rule for the same question:

if (material_index < 1)   // Material::Material

AT_material_index_from_material_number returns -1 for "not found" and 0 for the first table row — which today is User_Defined_Material = 0 (AT_DataMaterial.h:43,169, all-zero properties, ready = false). So < 1 happens to reject material 0 only because "User defined" is row 0 of AT_Material_Data; reorder that table upstream and this silently rejects a real material. The ID-based rule in the header is the correct one — have the constructor delegate to it:

if (id < 1 || AT_material_index_from_material_number(id) < 0)
  throw std::invalid_argument("Material not found: " + std::to_string(id));

Also: this is a breaking change that isn't called out. On master the check is material_index < 0, and index 0 is material 0, so pyamtrack.materials.Material(0) currently constructs and hands back an all-zero material. After this PR it raises ValueError: Material not found: 0. I think that is the right change — but it deserves a line in the PR description/release notes and a test (tests/test_materials.py currently only covers 9999).

H3. stopping.electron_range still has the old, contradictory material contract

The "single point of material validation" comment is aspirational: electron_range still routes through get_id(material, process_material), so the same user error gets two different answers depending on the module:

call	result
let_wilkens(1.0, 150.0, 0.01, 0)	ValueError: invalid material ID: 0
electron_range(120, 0)	returns nan — no error
let_wilkens(1.0, 150.0, 0.01, np.int64(1))	TypeError: material must be an int, a Material, ...
electron_range(120, np.int64(1))	RuntimeError: An error occurred while processing the Material object: ...

naming_convention.md presents material as one convention with one rule, so this divergence directly undercuts the doc. Migrating electron_range is a two-line change (validate_material_argument(material) + parse_material_argument(material) in place of get_id(...)) and would let process_material/get_id be deleted; if you'd rather not widen the PR, please soften the header comment and open a follow-up issue. Either way, the new helpers duplicate process_material + get_id rather than replacing them — the codebase now has two material paths where it had one.

H4. NumPy integer scalars are rejected, while NumPy integer arrays work

pm.let_wilkens(1.0, 150.0, 0.01, np.array([1, 2]))  # OK
pm.let_wilkens(1.0, 150.0, 0.01, np.int64(1))       # TypeError
pm.let_wilkens(1.0, 150.0, 0.01, True)              # OK -> silently material 1

np.int64 is not a PyLong, so nb::isinstance<nb::int_> is false and it is neither a list nor an ndarray — it falls through to the type_error. This is easy to hit (ids[0], arr.max(), anything indexed out of an int array). PyIndex_Check accepts exactly the integer-like objects and rejects floats:

} else if (PyIndex_Check(argument.ptr()) && !nb::isinstance<nb::bool_>(argument)) {
  id = nb::cast<long>(nb::steal(PyNumber_Index(argument.ptr())));
}

That also closes the material=True hole. Note the same fix is needed in parse_material_argument so the two stay in agreement, and it should return nb::int_ rather than the original object.

H5. let_wilkens returns NaN past the beam range, undocumented

pm.let_wilkens(17.0, 150.0)  # 24.88 keV/µm
pm.let_wilkens(18.0, 150.0)  # nan

Over np.linspace(0, 20, 2001) at 150 MeV with the default spread, 13.6 % of the returned values are NaN, silently. dose_bortfeld by contrast decays cleanly to 0.0 and stays finite everywhere I checked.

I pinned down the onset — it is not simply "past the range". In AT_LET_t/d_Wilkens_keV_um_single both nominator and denominator are clamped with if (x < 0) x = 0.0; and no zero guard follows, so once the Gaussian tail underflows the result is 0.0 / 0.0. Measured, the first NaN sits at exactly z − R = 5σ in every case, where R = α·E^p and σ = sqrt(σ_mono² + (σ_E·α·p·E^(p−1))²):

energy_spread_fraction	σ [cm]	first NaN (150 MeV, R = 15.693 cm)	(z−R)/σ
1e-6	0.1575	16.481 cm	5.00
0.001	0.1599	16.493 cm	5.00
0.01 (default)	0.3181	17.284 cm	5.00
0.05	1.3907	22.647 cm	5.00

So the usable depth window silently depends on energy_spread_fraction — a sharper-beam request shrinks it. That is worth stating in the Returns section rather than leaving users to discover a NaN inside plt.plot or an optimiser. Options, in increasing order of effort: document it; or return 0.0 past the range in the binding; or fix the 0/0 upstream in libamtrack and treat this as a known-issue note here. At minimum the example notebooks should keep their depth grid inside the finite window (examples/let_wilkens.ipynb currently goes to 20 cm at 150 MeV).

---

Medium

M1. export_values() leaks DOSE/TRACK into the module namespace

>>> [x for x in dir(pyamtrack.proton_models) if not x.startswith('_')]
['Averaging', 'DOSE', 'TRACK', 'dose_bortfeld', 'let_wilkens']

pyamtrack.proton_models.DOSE is now public API, undocumented in naming_convention.md (which lists only Averaging.DOSE/Averaging.TRACK). export_values() is the pybind11-era idiom for unscoped enums; for an enum class it just pollutes. Drop it. Also move the nb::enum_ registration above the two m.def calls — right now the type used by averaging is registered after the functions that accept it, which reads backwards.

M2. valid_material_ids is mutable, undocumented in the inventory, and untested

>>> pyamtrack.materials.valid_material_ids.add(999)   # succeeds, and persists process-wide

Three fixes: expose a frozenset (build the nb::set, then nb::steal(PyFrozenSet_New(s.ptr()))); add it to the pyamtrack.materials inventory in naming_convention.md (the id entry references it, but the inventory at §"Current Python API" doesn't list it); add a test. Also worth asking whether it earns its place at all next to the existing get_ids() — two public spellings of the same data is what the new doc is trying to prevent.

M3. Header declarations disagree with the definitions, and carry dead defaults

dose_bortfeld.h:26-28 declares z_cm, fluence_cm2, E_MeV, sigma_E_fraction; the definition uses depth_cm, ..., energy_MeV, energy_spread_fraction. Legal C++, but the header is what a reader reaches for, and its doc comment then documents names that don't exist in Python.

Worse, the declarations carry C++ default arguments — nb::int_(1), nb::float_(0.03), nb::cast(std::string("dose")) — that are (a) unreachable, since nb::arg(...) = ... in the bindings supplies the real defaults, and (b) a footgun: a default argument is evaluated at each call site, so these construct Python objects and require a live interpreter. Drop them and align the names with the definitions.

M4. Runtime error messages in dose_bortfeld.cpp use the internal names

Copilot flagged this for the docstring's Raises block; the same bug is in the messages a user actually sees:

ValueError: sigma_E_fraction must be in (0, 1), got: 0.000000
ValueError: z_cm must be >= 0, got: -1.000000
ValueError: E_MeV must be in [0.1, 10000.0], got: 0.000000

None of sigma_E_fraction, z_cm, E_MeV is a parameter of pyamtrack.proton_models.dose_bortfeld. let_wilkens.cpp gets this right — please match it.

M5. averaging is case-sensitive against its own enum

Averaging.DOSE is spelled uppercase, but averaging="DOSE" raises ValueError. A std::transform(..., ::tolower) before comparing removes a pointless trap.

M6. The docstring guarantees LET_d >= LET_t always; the implementation doesn't

Over np.linspace(0, 20, 2001) at 150 MeV, 358 of 2001 samples have LET_d marginally below LET_t, up to 1.07e-10 relative. That's an upstream numerical artifact at small spectral width, not a physics error — but the docstring states it as an absolute guarantee, and a test written from it fails. Either soften the wording ("up to numerical precision") or state a tolerance.

M7. validate_material_argument / parse_material_argument duplicate their dispatch

The two functions carry the same four-branch type dispatch, differing only in what they do at the leaf, and both are called on every request — so a NumPy material array is traversed three times (tolist() for validation, tolist() again for parsing, then the wrapper's own pass). One function that validates as it converts and returns the parsed object is smaller and strictly faster. Small nit inside it: nb::cast<long>(nb::cast<Material>(argument).id) resolves to the C++→Python overload only because .id on a temporary is an xvalue — nb::int_(nb::cast<Material>(argument).id) says what it means.

Relatedly, nb::type_error messages differ between the two functions ("an int, a Material, or a list / NumPy array of either" vs "an integer, Material, list, or NumPy array") for the same condition; which one a user sees depends on which function trips first.

M8. docs/naming_convention.md — good doc, a few claims don't match the code

• Multi-argument functions return np.ndarray for list input, not a list (verified). §"Scalar, list, and NumPy arguments" is ambiguous here, and the existing electron_range docstring is outright wrong ("a Python list for a list input") — this doc is the right occasion to fix that.
• "arrays must be C-contiguous" is unconditional, but only the cartesian path enforces it: pm.let_wilkens(np.arange(0, 10, dtype=float)[::2], 150.0) works fine elementwise.
• valid_material_ids is referenced under id but missing from the API inventory (M2).
• Nothing links to the file — docs/README.md links pre-commit.md and tests.md but not this. Worth adding, or the doc will be forgotten.
• Should mention that material accepts bool today (H4) — or better, stop accepting it.

---

Low / nits

• proton_models_bindings.cpp:116 fails clang-format (tab indent), i.e. pre-commit wasn't run on this branch. CI never checks it — .github/workflows/main.yml runs only build+pytest. Adding a pre-commit run --all-files job would have caught this and the CMakeLists.txt indentation automatically; worth doing separately.
• proton_models_bindings.cpp includes nanobind/stl/map.h, stl/vector.h and stl/string.h, and uses none of them — the file contains no std::map, std::vector or std::string at all.
• Material temp(ids[i]); m.attr(names[i].c_str()) = temp; — the temporary adds nothing. The loop is also bounded by ids.size() while indexing names[i]; they agree today (both skip table row 0 → 24 entries), but an explicit assert/size check documents the coupling for free.
• materials.cpp — the three explanatory comments deleted from to_name() are unrelated churn; please restore them (they explain non-obvious std::remove_if/std::transform lines).
• (long)mat_id — the codebase otherwise uses static_cast<>.
• dose_bortfeld docstring: energy_spread_fraction is missing optional and Default: 0.01 although the binding defaults it; the same entry mixes public and internal names ("sigma_E_MeV = energy_spread_fraction * E_MeV").
• let_wilkens docstring says "Calls AT_LET_d_Wilkens_keV_um" — the actual call is ..._single.
• dose_bortfeld.h:13 — "and optionally other args if you pass arrays" is informal for a public header.
• examples/dose_bortfeld.ipynb: np.linspace(0.0, 35.0, 800) puts >half the plot beyond a 150 MeV range (peak at 15.4 cm) — 0..20 shows the curve. The cells also use the old names as local variables (z_cm, E_MeV, sigma_E_fraction) and pass everything positionally, which undercuts the naming doc; depth_cm=..., energy_MeV=... would double as documentation. And no markdown narrative, unlike examples/example.ipynb.
• matplotlib is not in requirements-dev.txt — which is exactly why the committed notebook's only saved output is a ModuleNotFoundError. Please add it (or an examples extra) so pip install -r requirements-dev.txt gives you a runnable examples environment.
• 30 commits including several merge commits — fine if you squash, but the description then needs updating (B3).
• Pre-existing, but this PR widens the exposure: CMakeLists.txt:129 pins libamtrack to GIT_TAG master. The new bindings depend on AT_ProtonAnalyticalModels.h symbols, so an upstream rename now breaks builds of released tags retroactively. Pinning a tag/SHA is worth a follow-up issue.

---

Suggested order of work

1. Fix examples/let_wilkens.ipynb, or drop it from this PR (B3); update the PR description.
2. Add tests/test_proton_models.py (B2) — it will fail on B1 and H1, which is the point.
3. Fix the wrapper's exception swallowing + leak (B1).
4. Relax the spread bound to [0, 1) (H1); align the Material constructor with the shared rule and note the Material(0) behaviour change (H2).
5. Accept NumPy integer scalars, reject bool (H4); document NaN past the range (H5); drop export_values() (M1); rename the internal names in error messages (M4).
6. Then the doc/naming cleanups (M2, M3, M8) and nits.

Items H3 (unifying electron_range) and the pre-commit-in-CI / libamtrack-pinning suggestions are reasonable follow-up issues rather than blockers for this PR.

---

Generated by Claude Code

[Copilot]

Pull request overview

Copilot reviewed 17 out of 17 changed files in this pull request and generated 7 comments.