pyazr — the Python Interface

pyazr drives AZURE2 from Python in-process: the R-matrix engine is compiled into a pybind11 extension module (_azure2) and a session is a real C++ AZUREAPI object living in the interpreter. There are no subprocesses, no sockets and no port bookkeeping — the full R-matrix engine, the same code the GUI uses, becomes callable from a script.

This is the route to anything the GUI does not offer directly: custom minimizers, external samplers, parameter scans, systematic studies over many model variants, and publication figures built from the model rather than from exported files.

Installation

pip install -e .          # from the repository root
pip install -e ".[all]"   # plus what the examples need

Core dependencies are NumPy, mpmath and SciPy. mpmath is imported when pyazr.transform loads, so it is required rather than optional; the samplers and plotting libraries the examples use are the extras.

AZURE2 itself is C++ and is not installed by pip. The engine is the pybind11 module _azure2, built by CMake with USE_API=ON (the default) into the pyazr/ package directory; import pyazr loads it from there (or from the pip-installed package data).

First steps

from pyazr import azure2
import numpy as np

with azure2("13N.azr") as azr:
    best = np.asarray(azr.params_rwa, float)   # free parameters
    chi2 = np.sum(azr.calculate_chi2_rwa(best))
    print(chi2)

The context manager releases the engine — several MB of compound nucleus and data — rather than waiting for the garbage collector; without it, call close(). A closed session raises if used again.

Note

One session per azure2() object, and any number of them. Each azure2() owns an independent engine, so several can coexist in one interpreter — useful for sweeping model variants. They do not interfere: each enters its own directory for the duration of a call and leaves again, so your own working directory is never changed.

But drive them from one thread. The engine is not reentrant and keeps process-wide state. For parallelism give each process its own session (see Running in parallel).

The .azr stores its output/, checks/ and data paths relative to itself, which is what cwd= handles; it defaults to the .azr’s own directory, so paths you pass pyazr are resolved from wherever you are.

Two parameter conventions

AZURE2 exposes its parameters in two spaces, and mixing them silently produces wrong answers.

Reduced width amplitudes (*_rwa methods) are the natural fit space and the only one with analytic derivatives. Default to these.

azr.params_rwa               # the free parameter vector
azr.calculate_chi2_rwa(x)    # chi-squared per segment
azr.calculate_rwa(x)         # cross section per segment
azr.residual_jacobian(x)     # residuals + analytic Jacobian

Physical parameters (level energies in MeV, partial widths in eV) are what parameters.out reports and what the GUI shows.

azr.params                   # physical vector
azr.calculate(x)             # cross sections from physical parameters
azr.transform_rwa(x)         # rwa -> physical

Both vectors contain only the free parameters, in .azr order.

Warning

The gamma field in a .azr <levels> block is not a reduced width amplitude — it holds the physical value: Γ in eV for an open particle channel, an ANC for a closed one, Γ_γ in eV for a photon channel. The two differ by factors of 10² to 10⁷, and a file written with the wrong one loads without complaint and is wrong. Convert with transform_rwa.

Inspecting a model

Every parameter carries what it physically is, so you never have to parse a name like width_1_2:

print(azr.level_scheme)          # pairs -> J-groups -> levels -> channels
print(azr.parameters.table())    # every parameter with its metadata
print(azr.datasets.table())      # per segment: file, reaction, E range, norm error
print(azr.pairs)                 # masses, charges, spins, separation energies

for w in azr.parameters.widths.free:
    print(w.name, w.jpi, "L =", w.L, "S =", w.S, "pair", w.pair)

Filtered views (.free, .energies, .widths, .norms, .shifts) and lookups (.by_level(...), .by_name(...), .by_physical_level()) are all available. A LevelKey prints as 5/2-#2@6.588MeV; (jgroup, level) is its identity, since AZURE2 restarts level numbering inside every J-group.

Chi-squared and derivatives

chi2   = np.sum(azr.calculate_chi2_rwa(x))
val, g = azr.chi2_and_grad(x)          # analytic gradient
r, J   = azr.residual_jacobian(x)      # sum(r**2) == chi2

residual_jacobian computes the level-energy, reduced-width and normalization columns for roughly the cost of two forward evaluations, which makes Gauss-Newton and Levenberg-Marquardt fits cheap. It also yields the per-segment χ² the API does not expose directly, by slicing the residuals by segment length.

Energy-shift columns are finite-differenced instead, at two extra residual evaluations each. A shift translates the energy axis of a whole segment, so what is wanted is \(\partial\,\mathrm{model}/\partial E\), and AZURE2 applies a shift by rebuilding every energy-dependent quantity of the affected points — level matrix, penetrabilities, Coulomb and hard-sphere phases, external-capture amplitudes, and the target-integration sub-point grid. A fit with many free shifts is dominated by their columns.

Warning

These return the data χ² only. AZURE2’s own fit objective additionally penalises every free normalization and energy shift. Minimise the residuals alone and the normalizations will drift to absorb every discrepancy, reaching a χ² that AZURE2 would never find. Append the penalty rows — see How the Chi-Squared Is Built and pyazr/examples/ for the recipe.

Editing the model

AzrModel parses the <levels> block and re-emits everything else verbatim, so the original file is never modified:

from pyazr import AzrModel

mdl = AzrModel.from_file("7Be.azr")
mdl.remove_level(jpi="1/2+", energy=20)
mdl.add_level(J=1.5, parity=+1, energy=8.6,
              channels=[dict(pair=1, L=2, S=0.5, gamma=1000.0, fixed=False)])
path = mdl.write("_variant.azr")

Because AZURE2 reads its model from the file, an edited scheme takes effect by opening a fresh azure2() session on the new file.

One rule catches people out: all levels of one J^π share a channel set. Adding a level to an existing J^π clones that group’s channel structure; only a brand-new J^π lets your channel list define the group.

Also available: set_channel_radius, set_segment_norm, set_segment_active, set_segment_datafile, add_data_segment, remove_data_segments, clear_data_segments, set_extrapolations and apply_fit, which writes a fit result back into a .azr you can reopen in the GUI.

Warning

Adding or removing data segments changes which energies AZURE2 evaluates. The external-capture integrals are cached in output/intEC.dat / output/intEC.extrap keyed on the grid, not on the segment selection, and AZURE2 silently reuses a stale file. After any data edit, delete those caches (or give the edited model its own output directory) so the integrals are recomputed; azure2.recalculate_external_capture() forces it inside a live session. See pyazr/examples/edit_model.py.

Decomposing a cross section

At a fitted point, without refitting, each number is a component’s raw contribution:

azr.extrap_mode()
full = azr.calculate_rwa(best)
off  = azr.calculate_rwa(azr.without_level(best, jpi="5/2-", energy=6.588))
only = azr.calculate_rwa(azr.only_level(best,   jpi="5/2-", energy=6.588))

full - off is everything the level does, resonance plus interference; only minus the non-resonant background is the bare resonance. The difference between them is pure interference, which is block-diagonal in J^π — only same-J^π levels interfere, and that is a useful check on any decomposition.

Fetching data from EXFOR and NDS

pyazr.nds wraps the IAEA web services for experimental nuclear data and converts the results straight into AZURE2 form (network required). EXFOR holds the measured reaction data; LiveChart/ENSDF the evaluated level schemes.

from pyazr import nds

# search EXFOR: target, reaction (proj,exit) and quantity filters
hits = nds.search_exfor(target="C-13", reaction="p,g", quantity="SIG")
# -> [ExforDataset O2599004 6-C-13(P,G)7-N-14,,SIG,,SFC n=31 ...]

# fetch one dataset, convert to an AZURE2 data file + segment kwargs
data = nds.fetch_exfor("O2599004")                 # S-factor, B*KEV
kw = data.to_azr("run/data", entrance=1, exit=2,
                 observable="total-capture")       # lab E, barns
AzrModel.from_file("13N.azr").add_data_segment(**kw).write("13N_new.azr")

# level scheme of the compound nucleus, from ENSDF
for lv in nds.fetch_levels("14n"):
    print(lv.energy_mev, lv.jp, lv.half_life)

# the paper behind the dataset, resolved to a DOI
ref = nds.reference("O2599004")
doi = nds.resolve_doi(ref)   # -> 10.1103/physrevlett.131.162701

Frames and units are handled automatically: EN-CM energies are converted to lab, B*KEV S-factors to barns via the Sommerfeld factor, NB/SR differentials to b/sr. Pick observable to match the data’s angle frame (EXFOR ANG-CM → "differential-cm" or "analyzing-power"; lab ANG → "differential"). See the nds-explorer skill for the full quantity/unit reference and failure modes.

One case cannot be decided automatically. Ratio-to-Rutherford data (,,RTH) is dimensionless with an angle column, which is exactly what an analyzing power looks like — and x4get drops the quantity suffix that distinguishes them, so only the search result knows. to_azr warns and passes such values through; give it rutherford=True to multiply by the Coulomb cross section, or rutherford=False once you have checked.

The module needs nothing beyond NumPy and the standard library, and doubles as a command-line tool:

python -m pyazr.nds search --target C-13 --reaction p,g --quantity SIG
python -m pyazr.nds download O2599004 -o data/skowronski.dat
python -m pyazr.nds reference O2599004

The Qt setup utility has its own EXFOR dialog, backed by gui/src/ExforData.cpp. The two are independent implementations of the same Web-API, so a parsing rule learned by either belongs in both.

Angular distributions

AZURE2 writes an angular distribution as

\[W(\theta) = \sum_k a_k P_k(\cos\theta)\]

and computes the \(a_k\) only for segments declared as angular distributions – observable="angular-distribution" with an order. Everything else yields empty arrays.

For the grids a model already declares, on a live instance:

dists = azr.calculate_angular_dists_rwa(x)
# one entry per segment; each a list with one array of coefficients per point

For arbitrary energies, which is usually what is wanted:

from pyazr import angular_distribution

e_cm, coeffs = angular_distribution("model.azr", [0.05, 0.1, 0.2],
                                    entrance=1, exit=2, order=4)

coeffs[i, k] is \(a_k\) at energy i. Input energies are lab by default (lab=False for centre-of-mass); the returned energies are always centre-of-mass. Energies AZURE2 could not evaluate come back as NaN rather than being dropped, so the rows always line up with the input.

Each call writes a temporary model requesting exactly those energies and opens one azure2() session on it, so pass every energy in a single call rather than looping.

Reading the result: \(a_0\) is the normalisation, so \(a_0 = 1\) with everything else zero means isotropic. That is the correct answer for a resonance formed in an s-wave – the compound nucleus has no preferred direction – and is what ³H+d gives at low energy. Anisotropy shows as the higher orders departing from zero.

Worked example: pyazr/examples/angular_distribution.py.

Analyzing power

The vector analyzing power \(A_y\) is observable code 7. Declare a segment with observable="analyzing-power" and it is reported in place of the cross section, so \(\chi^2\), fitting and plotting need no special handling:

ay = azr.calculate_analyzing_power_rwa(azr.params_rwa)

Data files carry E_lab theta_cm A_y dA_y – angles are centre-of-mass, unlike an ordinary differential segment. Leave vary_norm off: a normalization factor is meaningless for a ratio.

One trap is worth knowing about: comparing against thin-target data requires segments with no target integration, because \(A_y\) averaged over a thick target is weighted by the cross section and Rutherford scattering drives that average towards zero.

residual_jacobian and chi2_and_grad do differentiate \(A_y\) exactly. The one exception is an analyzing-power point that also carries target integration: that is a ratio of two integrals, is not differentiated analytically, and makes the analytic Jacobian unavailable for the whole fit rather than returning something approximate.

See Polarization Observables in R-Matrix Theory for the formalism, Computing the Analyzing Power in AZURE2 for the implementation, and tests/13N (segments 11–16) for a worked comparison against measured data.

Cleaning up

There are no stray processes to reap: a session runs in-process, so when the interpreter exits, the engine goes with it. No orphans, no pyazr.cleanup, no port collisions.

There is still memory. A session holds its compound nucleus and data — several MB for a small model, more for a real one — so in a loop over model variants use the context manager (or close()) rather than leaving it to the garbage collector, which may otherwise hold two models at once.

Dimensionless widths

t = azr.dimensionless_widths(best)
print(t.particles.table())         # theta^2 per particle channel
print(t.photons.nonzero.table())   # Weisskopf units per gamma channel
[c for c in t.particles if c.theta2 and c.theta2 > 1]   # unphysical

θ² > 1 exceeds the Wigner limit and is unphysical; γ-ray strengths in Weisskopf units are the corresponding sanity check for capture.

The external region, and the caches

The quantities that describe everything outside the channel radius can be asked for directly. They are what the penetrabilities in the level matrix are built from, what sets the hard-sphere phase, and what the external-capture integrals integrate.

c = azr.coulomb_functions(pair=1, energies=E, L=0)   # radius=0 -> channel radius
c["F"], c["G"], c["P"], c["S"], c["delta_hs"]

paths = azr.ec_integrals(pair=1, energies=E)         # one entry per EC pathway
paths[0]["li"], paths[0]["lf"], paths[0]["radiation"], paths[0]["value"]

azr.cache_stats()   # queries, hits, hit_rate, entries, keys, disabled_keys

The Coulomb functions follow the run’s own configuration, so the same call returns the accurate routine’s values, GSL’s (--gsl-coul), or the Numerov solution through a nuclear potential (the hybrid model of the <potential> block). Comparing them is how one sees what those options do to the external region.

External-capture integrals are the most expensive part of a capture calculation, which is why the Coulomb functions they need are memoized. cache_stats makes that visible: asking for the same integrals twice on 16O(p,gamma)17F takes 14.9 s and then 0.44 s, with the hit rate rising from 82% to 91%. disabled_keys counts the memos that have given up because too few of their entries were being asked for twice — which is what a varying energy shift produces, since it moves every point energy at every iteration.

Examples

Worked scripts ship in pyazr/examples/:

Script	What it shows
`print_scheme.py`	Reading a model’s level scheme and dataset provenance.
`edit_scheme.py`	Adding and removing levels, writing a new `.azr`.
`edit_model.py`	The full edit loop: add/remove resonances and data segments, recalculate the external-capture integrals when the data change, and save the edited file.
`exfor_fetch.py`	Pulling data from the IAEA EXFOR/NDS web services (cross sections, S-factors, differentials, analyzing powers, levels) and dropping it into a model as new segments — with the DOI of the paper behind the dataset.
`deactivate_level.py`	Switching a resonance off without removing it.
`transform_widths.py`	Reduced width amplitudes to physical partial widths.
`dimensionless_widths.py`	θ² and Weisskopf units for a whole fit.
`coulomb_functions.py`	Coulomb functions, penetrability and hard-sphere phase over an energy grid.
`ec_integrals.py`	External-capture integrals per pathway, and what caching them buys.
`save_fit_to_azr.py`	Writing a fit back into a `.azr`, verifying it round-trips.
`uncertainty_band.py`	Cross-section uncertainty bands from a saved fit covariance.
`sensitivities.py`	Which parameters a dataset constrains, as d ln σ / d ln p; and the analytic sensitivities checked against finite differences.
`fit_emcee.py`	MCMC sampling with `emcee`, one in-process engine per pool worker.
`fit_zeus.py`	The same with `zeus`.
`angular_distribution.py`	Legendre coefficients at chosen energies, with an optional plot.
`per_dataset_chi2.py`	Slicing the residual vector to get the χ² of each experiment.
`decompose_cross_section.py`	Separating resonance, interference and background contributions at the fitted parameters.
`sfactor_extrapolation.py`	An S factor below the measured range, and S(0).
`channel_radius_scan.py`	χ² against the channel radius, one in-process session per radius.
`reaction_rate.py`	N_A⟨σv⟩ by integrating the extrapolated cross section over a Maxwell–Boltzmann distribution.

Running in parallel

There are no AZURE2 instances to pool — every azure2() object is already an independent engine in-process. But the engine is not thread-safe, so the unit of parallelism is the process: give each worker its own session by constructing it at module scope in the worker module. That works whichever way the pool starts its workers — under spawn each re-imports the module and builds its own engine, under fork each inherits a copy — and each then evaluates one χ² per walker:

# worker module
azr = azure2("13N.azr")

def log_prob(theta):
    return -0.5 * np.sum(azr.calculate_chi2_rwa(theta)) + prior(theta)

See pyazr/examples/fit_emcee.py and fit_zeus.py for the full recipe.

The same applies inside pyazr itself: sensitivities, uncertainty_bands and extrapolation_bands take nprocs=N, which spreads the finite-difference columns over N worker processes, each with its own engine and its own output directory (seeded from yours, so the external-capture integrals are reused rather than rebuilt).

It is not free, and it is not always a win. Every worker pays one model initialisation before it evaluates a single column, so the columns have to be worth more than that startup — which in practice means hundreds of free R-matrix parameters. On the 14-parameter tests/13N model it loses at any nprocs, both in data mode (0.8 s serial) and on an 1845-point extrapolation grid (53 s serial, 56 s over four workers). It is also irrelevant to the default method="analytic", which is one call whatever the parameter count. Reach for it only when the analytic path is unavailable and the model is large, and measure rather than assume.

Data mode and extrapolation mode

data_mode() (the default) evaluates the <segmentsData> segments — this is what χ² uses. extrap_mode() switches to the <segmentsTest> grids for predictions on arbitrary energies and angles. Both re-initialise the engine, so switch sparingly and batch the work.

Segment indexing differs between them, and getting it wrong misaligns results silently: in data mode segment i is azr.datasets[i]; in extrapolation mode only active test segments are returned, so segment i is azr.extrapolations.active[i].

Note

Delete output/intEC.extrap whenever the <segmentsTest> grid changes. AZURE2 caches external-capture integrals there and silently reuses them on a different grid, which corrupts capture cross sections. It is safe to delete; it only costs time to rebuild.

Frames and units

Input is lab frame — the energies and angles in a .azr and in data files. All output and every API result is centre-of-mass. This includes extrapolation grids: the energies you set are lab, the energies you get back are c.m. Convert with E_lab = E_cm · (m_beam + m_target) / m_target.

Plot different reaction channels against excitation energy (calculate_excitation_energy); it is the only axis shared by all entrance channels.