pyextal.jaxMultipole.intensity¶
Option 2: multipole parameters refined directly against LARBED intensities.
The joint refinement (jaxJoint) refines one table of shared structure factors
– G['shared'], one {xfre, xfim, afre, afim} per canonical family – that
every dataset reads through its own sign/phase links. Here that table is no
longer free: its elastic part is computed from the multipole model,
xfre + i xfim = U_g(p), p = p0 + M (s * x),
for every shared family, in float32, inside the differentiable forward model. The free parameters are
xthe multipole parameters chosen by a ParameterMap (populations, kappa, P_lm, … – whatever it frees), scaled by
s(typically their Option 1 esd) so that one unit is a comparable step for Adam;absorptionthe absorptive U’_g of each shared family, still free (
afre;afimas the joint run’simagoption says), as in the joint refinement;datasetseach dataset’s thickness, scale, background and probe, unchanged.
Everything else is jaxJoint’s: the same per-dataset losses, links, chunking and pooled weighting, and between Adam segments the same Levenberg-Marquardt geometry refinement of each dataset.
Scope (step 2A)¶
Only the shared families come from the model. The other Bloch-matrix elements
– the differences between beams that no dataset refines – keep the engine’s
tabulated IAM (Cromer-Mann, the .dat Debye-Waller factors), as they do in
the joint refinement; so do the beam selection (beamcheck) and U_0.
Warm start¶
warm_start copies a finished joint run into the tree: each dataset’s refined
geometry (from joint_<stamp>.json), its thickness, scale, background and
probe (from its run_NNN in the dataset file) and the shared absorption. The
multipole start comes from the Option 1 fit. The first loss is then the joint
run’s, except that the elastic U_g are now the model’s.
Uncertainties¶
curvature_esd evaluates the Hessian of the pooled loss with respect to x.
The pooled loss is the mean chi^2 per data point, so chi^2_total = N L and
cov(x) = 2 (N H)^-1; scaled by L (the reduced chi^2, well above 1 for these
data) it is the conventional esd. Correlated pixels (the probe convolution, the
geometry) make even that optimistic.
Classes¶
Options of refine_intensities beyond the jaxJoint ones. |
|
Maps the multipole parameter tree onto jaxJoint's |
|
The multipole model on every coupling of the Bloch-wave matrix (Option 2B). |
Functions¶
|
Copies a finished joint run's geometry, scalars and absorption into |
|
Damped Newton on |
|
Curvature and dataset-jackknife esd of |
|
Loss cost of each multipole U_g in the joint run's own solution. |
|
Datasets, the warm-started joint tree and the multipole bridge. |
|
One dataset's simulated and measured rocking curves for the tree |
|
Per reflection: R = sum|sim - exp| / sum|exp| and rms of (sim - exp) / max(exp), in %. |
|
Experiment | reference simulation | new simulation | new - reference, per reflection. |
|
Simulated rocking curves of the multipole solution |
|
Option 2 in stages: each stage refines its own parameters, the rest frozen. |
|
Option 2: the multipole model refined against the joint LARBED intensities. |
Module Contents¶
- class pyextal.jaxMultipole.intensity.IntensitySettings¶
Options of refine_intensities beyond the jaxJoint ones.
- Variables:
mode (str) –
'newton'(default) – damped Newton on the multipole parameters with everything else held, using the exact Hessian of the pooled loss (a few parameters: a few gradient evaluations per iteration), optionally followed byrest_stepsAdam steps on the absorption and per-dataset scalars, repeatedcyclestimes;'adam'– Adam on everything at once forstepssteps, with jaxGeo’s scaling and warm-up (seesteps).newton_lambda (newton_iters, newton_tol,) – Newton iterations, the step size (in units of x) below which it stops, and the initial damping.
hessian_every (int) – Accepted Newton steps between exact Hessians (
len(x)Hessian-vector products, ~2-3 gradients each).rest_steps (int) – Adam steps on absorption + scalars after each Newton stage.
cycles (int) – Newton / Adam alternations.
family_cost (bool) – Before fitting, report the loss cost of each model U_g in the joint solution (family_cost).
steps (int) – Adam steps (
mode='adam'): the multipole parameters, the absorption U’_g and the per-dataset scalars together, scaled as jaxGeo scales them – onelr; the scalars normalised to their windows, the absorption in A^-2, and each multipole parameter’s step scaled (MultipoleBridge.ug_units) so that one step moves the model U_g by at mostlr_xA^-2, as one step moves a free U_g in jaxGeo.lr (float) – Adam step (jaxGeo’s
lr).lr_x (float) – The U_g step of the multipole parameters, A^-2 (
lrwhen None).warmup (int) – Steps at the start with the multipole parameters and the absorption frozen, then a fresh Adam state – jaxGeo’s
warmup; the joint run’s value when None.polish (bool) – After Adam, damped Newton on the multipole parameters at the final absorption and scalars.
lr – Adam step for absorption and per-dataset scalars (as jaxJoint).
geometry_every (int) – Geometry LM every this many steps (0: never).
esd (bool) – Evaluate the curvature esd at the end.
save_runs (bool) – Simulate and save a
run_NNNin every dataset file.out (Path) – Where the summary JSON goes.
log_every (int) – Print every this many steps.
- mode: str = 'newton'¶
- newton_iters: int = 12¶
- newton_tol: float = 0.005¶
- newton_lambda: float = 0.001¶
- hessian_every: int = 3¶
- rest_steps: int = 0¶
- cycles: int = 1¶
- steps: int = 200¶
- lr: float = 0.0005¶
- lr_x: float = None¶
- warmup: int = None¶
- polish: bool = False¶
- geometry_every: int = 0¶
- esd: bool = True¶
- save_runs: bool = False¶
- out: pathlib.Path = None¶
- log_every: int = 5¶
- family_cost: bool = False¶
- class pyextal.jaxMultipole.intensity.MultipoleBridge(model, pmap, keys, scale=None)¶
Maps the multipole parameter tree onto jaxJoint’s
G.- Parameters:
model (MultipoleModel) – Its crystal must be the one the datasets describe (same cell, setting and voltage).
pmap (ParameterMap) – Which multipole parameters are free.
keys (list) – The shared canonical keys, in
G['shared']order.scale (array, optional) – Units of
x(one per free parameter); 1 by default.
- keys¶
- scale¶
- geo¶
- _p0¶
- _Ms¶
- p(x)¶
Physical parameters for scaled
x(float32).
- ug(x)¶
Model U_g of the shared keys, complex64.
- G(theta)¶
jaxJoint’s
Gfor the treetheta.
- theta_from(G, x0=None)¶
The tree for a joint
G: its absorption and scalars,x0(zeros).
- ug_units(x=None)¶
Per free parameter: the step in
xthat moves the shared U_g by at most 1 A^-2.1 / max_g |dU_g / dx_i|atx(zeros by default). Multiplying an Adam update bylrtimes this makeslrmean what it means for a free U_g in jaxGeo / jaxJoint: the largest U_g change per step, in A^-2.
- physical(x)¶
- class pyextal.jaxMultipole.intensity.FullMatrix(model, datasets, iam_model=None, log=print)¶
The multipole model on every coupling of the Bloch-wave matrix (Option 2B).
MultipoleBridge drives only the shared, refined structure factors; every other U_(g-h) in a dataset’s matrix – the beams that are coupled but not refined, upper-layer (HOLZ) couplings such as (111), (211) included, and the weak beams of the Bethe sums – stays at the engine’s static IAM (Cromer-Mann with the
.datDebye-Waller factors). This class adds, at every one of those positions, the changeU_model(g) - U_IAM(g)as the engine’sugh_delta(jaxEngine.model.construct_ugh), so the whole matrix follows the multipole parameters and the refined ADPs.Set-up, per dataset: the distinct difference vectors
g_i - g_jof its beam list, each mapped to its canonical family with the symmetry link (U(g) = c U(key)orc U(key)*, symmetry.link_phase); the positions not covered by a refined reflection (metadata['mask']). The model is evaluated once per call for the union of families of all datasets.The absorptive part of these couplings stays the engine’s Bird-King (
.datB); the refined reflections keep their own route (model elastic part, refined U’_g).- Parameters:
model (MultipoleModel) – The multipole model (its crystal = the datasets’).
datasets (list) – jaxJoint datasets (after prepare).
iam_model (MultipoleModel, optional) – The IAM the static matrices were built from; by default the model’s crystal in
'cm'mode, which reproduces the engine’sdefaultsf. check compares it with the static matrices.
- model¶
- per¶
- keys = []¶
- geo¶
- iam_model¶
- U_iam¶
- dU(p)¶
U_model - U_IAMof every family, complex (pphysical, float32).
- delta(name, p=None, dU=None)¶
The engine’s
ugh_deltaof datasetname: per layer(dr, di).
- iam_matrix(name)¶
The IAM elastic matrix this class assumes the static one holds, (nb, nb) complex.
- check(datasets, log=print)¶
Compares the assumed IAM with each dataset’s static matrix (real part, free positions).
- Returns:
name -> (max |diff|, max |static|)in A^-2.- Return type:
dict
- pyextal.jaxMultipole.intensity.warm_start(datasets, G, joint_json, engine)¶
Copies a finished joint run’s geometry, scalars and absorption into
G.- Parameters:
datasets (list) – jaxJoint Dataset objects (after prepare).
G (dict) – The joint tree from prepare; modified in place.
joint_json (path) –
joint_<stamp>.jsonof the run to start from.engine (EngineSwitch) – For the geometry rebuilds.
- Returns:
What was copied, per dataset.
- Return type:
dict
- pyextal.jaxMultipole.intensity.newton_x(theta, work, iters=12, tol=0.005, lam=0.001, hessian_every=3, log=print, label='', bounds=None)¶
Damped Newton on
theta['x']with everything else held.Solves
(H + lam diag|H|) dx = -g; a step is accepted when the pooled loss falls (thenlamshrinks) and otherwiselamgrows. The exact Hessian (len(x)Hessian-vector products per chunk) is recomputed everyhessian_everyaccepted steps, or when a step with a stale one fails; between, only the gradient is. Stops when the accepted step is belowtol(units of x) or afteritersiterations.bounds((lower, upper)in units of x) clips every trial step, so a radial parameter cannot leave its physical range (e.g. kappa’ < 0.3).- Returns:
(theta, history[(loss, x, per-dataset)], H).- Return type:
tuple
- pyextal.jaxMultipole.intensity.curvature_esd(theta, work, n_points)¶
Curvature and dataset-jackknife esd of
xat the optimum.Two estimates:
formal:
cov = 2 (N H)^-1from the HessianHof the pooled loss (a mean overNpixels). It treats every pixel as an independent measurement; the probe convolution, the tiling and the model error make neighbouring pixels strongly correlated, so it is far too small;jackknife over datasets: each dataset’s own Newton step from the optimum gives the delete-one estimates
x_(-d) = x - (H - w_d H_d)^-1 (g - w_d g_d); their spread is the esd that the scan-to-scan reproducibility supports (the analogue of the per-scan scatter used for the Option 1 sigma). Needs >= 3 datasets.
- Returns:
hessian,loss,cov,esd(formal),esd_scaled(formal times sqrt L),corr,jackknife(esd or None),per_dataset_x(each dataset’s own optimum, one Newton step), all in the units ofx.- Return type:
dict
- pyextal.jaxMultipole.intensity.family_cost(datasets, G, bridge, x, log=print)¶
Loss cost of each multipole U_g in the joint run’s own solution.
For every shared family, the joint run’s refined
U_gis replaced by the model’s (all else held at the joint solution) and the pooled loss re-evaluated: the increase is what the intensities say against the model’s value of that family. Costs one loss evaluation per family and dataset.- Returns:
(key, U_joint, U_model, cost)sorted by cost, largest first.- Return type:
list
- pyextal.jaxMultipole.intensity.setup_state(js, model_builder, pmap_builder, joint_json=None, scale=None, log=print)¶
Datasets, the warm-started joint tree and the multipole bridge.
- Returns:
model,pmap,datasets,G(the joint run’s solution, its own refined U_g),table,engine,bridge,theta(x = 0: the pmap start),keys,js(resolved).- Return type:
dict
- pyextal.jaxMultipole.intensity.simulate_dataset(ds, G, ugh_delta=None)¶
One dataset’s simulated and measured rocking curves for the tree
G.As jaxJoint.refine._finish draws them (over-subtraction applied, the experiment clipped to match), but nothing is written to the dataset’s HDF5.
- Returns:
(sim, exp, probe)– arrays (rows, cols, n) and the Voigt probe dict for save_comparison_png (or None).- Return type:
tuple
- pyextal.jaxMultipole.intensity.reflection_stats(sim, exp)¶
Per reflection: R = sum|sim - exp| / sum|exp| and rms of (sim - exp) / max(exp), in %.
- pyextal.jaxMultipole.intensity.save_three_way_png(exp, sim_ref, sim_new, labels, path, names=('joint (free U_g)', 'multipole'), stats=None)¶
Experiment | reference simulation | new simulation | new - reference, per reflection.
The last column, in % of the experiment’s peak on a shared colour scale per reflection, shows where the multipole U_g change the intensities relative to the joint run’s own refined U_g, with every other parameter identical.
- pyextal.jaxMultipole.intensity.draw_intensities(state, x, outdir, tag, names=None, joint=True, family_of=None, log=print)¶
Simulated rocking curves of the multipole solution
xfor each dataset.Writes, per dataset,
<name>_<tag>_comparison.png(the jaxGeo layout: experiment, simulation, residual in counts and in %, line profile) and, withjoint,<name>_<tag>_vs_joint.png(experiment, the joint run’s own simulation, the multipole one, and their difference). Nothing is written to the datasets’ HDF5 files. Withfamily_of(hkl -> family key) the figures show one reflection per family – its member with the most counts; the statistics always cover every reflection.- Returns:
per dataset
dict(name, labels, R_joint, rms_joint, R, rms).- Return type:
list
- pyextal.jaxMultipole.intensity.refine_sequence(js, model_builder, stages, start, joint_json=None, settings=None, theta_init=None, final_refine=None, scale_of=None, out=None, full_matrix=False, log=print)¶
Option 2 in stages: each stage refines its own parameters, the rest frozen.
The datasets are set up once; every stage builds a parameter map from the current physical values (so frozen parameters stay exactly where the previous stages left them) and runs damped Newton (newton_x) on the joint intensities, with absorption and per-dataset scalars held.
- Parameters:
stages (list) –
(name, refine)withrefinea build_parameter_map refine dict; a stage whose map has no free parameter is logged and skipped (e.g. l = 1, 3 on centrosymmetric sites).start (dict) – Physical start values (names of
model.layout).final_refine (dict, optional) – Everything free at once, evaluated at the end point without steps: the marginal jackknife esds (the per-stage ones are conditional on the frozen parameters, hence too small).
scale_of (callable, optional) –
name -> unit of x(Newton is invariant to it except through the damping and the stopping step).out (path) – Summary JSON (the refine_intensities keys for the end point, plus
stages); the tree goes next to it as.theta.pkl.full_matrix (bool) – Drive every coupling of the Bloch-wave matrix by the model (FullMatrix), not only the refined structure factors.
- Returns:
stages(per stage: names, values, loss, per-dataset losses, time),values(final physical),esd(iffinal_refine).- Return type:
dict
- pyextal.jaxMultipole.intensity.refine_intensities(js, model_builder, pmap_builder, joint_json=None, settings=None, scale=None, theta_init=None, full_matrix=False, log=print)¶
Option 2: the multipole model refined against the joint LARBED intensities.
- Parameters:
js (JointSettings) – The datasets and common options, as for jaxJoint.
model_builder (callable) –
() -> MultipoleModel; called before the datasets are set up (building a Crystal re-parses the engine).pmap_builder (callable) –
model -> ParameterMap(free parameters and the physical start).joint_json (path, optional) – Warm start from this finished joint run.
settings (IntensitySettings) – Steps, learning rates, geometry, output.
scale (array, optional) – Units of
x; e.g. the Option 1 esd.theta_init (dict, optional) – A previous run’s final tree (its
.theta.pkl): its absorption and per-dataset scalars replace the joint run’s; itsxis ignored (the physical start comes frompmap_builder).
- Returns:
theta,bridge,datasets,history,per_dataset,esd(if requested),summary(path).- Return type:
dict