diff --git a/docs/SCHEMA_HISTORY.md b/docs/SCHEMA_HISTORY.md index 0e31e2b..7667bf1 100644 --- a/docs/SCHEMA_HISTORY.md +++ b/docs/SCHEMA_HISTORY.md @@ -192,9 +192,65 @@ engine schemas adopt it in re/04–06. ## gsasii engine schema -### 1.0.0 — introduced in PowderLine 0.2.0 +### 1.0.0 — introduced in PowderLine 0.2.0 (in development) -Entries added by re/04 (`gsasii.rietveld`, `gsasii.spf`). +Defined in `src/powderline/gateways/gsasii/schema.py` (engine-free), on core +1.0.0. Accepted versions are declared: engine schema `==1.0.0`, requiring +core `==1.0.0`. A 0.26.0 `GSASII_*` recipe is converted with +`scripts/convert_recipe_026.py`, which reports every change. + +- **Workflows**: `gsasii.rietveld` (`phases`, no `single_peaks`) and + `gsasii.spf` (`single_peaks`, no `phases`; `refinement_controls.single_peak_fitting_mode` + required). +- **Parameters**: `[value, refine_flag]` only. GSAS-II applies parameter + limits only after the fit, so gsasii 1.0.0 exposes no bounds; a 4-element + list is an error. +- **Phase block**: core's `space_group`, `unit_cell`, `atoms`, plus `scale` + (required, >= 0) and `peak_broadening` (isotropic size in µm and microstrain, + each with `LG_eta`). The phase name is its key in `phases`. + - Space group: GSAS-II cannot use origin choice 1 or 9 other settings (7 it + does not read; `P 21 n m` and `A b a m` it reads as another setting), so + they are errors. + - Cell refinement: GSAS-II refines an oblique cell's parameters together + (monoclinic {a, c, β} for b-unique, {b, c, α} a-unique, {a, b, γ} + c-unique; all six for triclinic and rhombohedral axes). Refining part of + such a group is an error. + - **Size/strain left out**: 10 µm / 0 microstrain, isotropic, fixed, with a + structured warning, instead of GSAS-II's own 1 µm / 1000 microstrain. A + stated crystallite size must be > 0. +- **Instrument** (single source; GSAS-II's instrument parameters are built from + it): `description` (optional), `radiation {type: "PXC", wavelength}`, + `geometry {bank, azimuth}` (plain values), `corrections {zero_shift, + polarization, axial_divergence}`, `broadening {U, V, W, X, Y, Z}`. GSAS-II + units: U, V, W centideg²; X, Y, Z and zero centideg; SH/L dimensionless. + SH/L below 0.002 is an error (GSAS-II computes with 0.002 anyway). The + wavelength must be > 0. Kα doublets, TOF and neutron are not supported. +- **Background**: `chebyshev` (left out: one fixed term 0.0, no background, + instead of GSAS-II's constant 1.0) and optional `single_peaks` (positions, + intensities, `pv_gaussian_sigma_sq` in centideg², `pv_lorentzian_gamma` in + centideg; all four lists the same length; each position inside the fit + window). Values below GSAS-II's silent floors (intensity 0.1, σ² 0.01, γ 0.1) + are errors. Peak List peaks (`gsasii.spf`): σ² and γ at least 0.001, each + position inside the fit window. +- **Single peak fitting widths**: with `use_instrument_profile: true` the peak + widths come from the instrument profile (U..Z may be refined; a peak's own + σ²/γ may not); with `false` each peak has its own widths (they may be + refined; U..Z may not). GSAS-II silently ignores a refine flag on the side + not in use, so such a flag is an error. +- **Simulation** (`refinement_cycles` 1): every refine flag must be false. +- **At run time**, after phase setup: an atom GSAS-II reads with the wrong + multiplicity (some 2-fold sites in R32, R-3m, R-3c), or whose site it cannot + name while a coordinate or Uaniso is refined, is an error naming an + equivalent position to state instead. +- **Results**: `rwp`, `r_exp`, `gof`, `chi2_red` from PowderLine's uniform fit + statistics over the stated fit window (`None` where undefined, e.g. no + degrees of freedom); `simulation_mode` (a simulation gets the same + statistics; against placeholder data they mean nothing); GSAS-II's own + values, `parameters_requested` and `parameters_varied` under + `engine_details`; structured `warnings` (defaults applied, GSAS-II's + refinement messages, fewer parameters varied than requested). A refinement + that diverges (a non-finite calculated pattern) is a failure with a clear + message; GSAS-II itself does not flag it. ## topas engine schema diff --git a/scripts/convert_recipe_026.py b/scripts/convert_recipe_026.py new file mode 100644 index 0000000..026865c --- /dev/null +++ b/scripts/convert_recipe_026.py @@ -0,0 +1,524 @@ +#!/usr/bin/env python3 +"""Convert a 0.26.0 recipe (``GSASII_Rietveld`` / ``GSASII_SPF``) to gsasii 1.0.0. + +The converted recipe states explicitly what v0.1.1 effectively refined (master +plan A95, A101, re/04 subplan §3): structure and parameterization merge into one +phase block, every refinable becomes ``[value, refine_flag]``, refine flags are +expanded to the groups GSAS-II refined together, and GSAS-II's own silent +defaults are written out. Every change is reported. Nothing is fixed silently: +a value the gsasii 1.0.0 schema rejects stops the conversion with every +problem listed (e.g. a background peak below GSAS-II's floor, an ambiguous +special position, a name outside the A112 rule). + +Engine-free (no GSAS-II import). Used by the re/04 equivalence runs and the +re/07 example migration; the user-facing migration guide is written in re/10. + +Usage:: + + python scripts/convert_recipe_026.py IN.json [-o OUT.json] + +API: :func:`convert` returns ``(recipe, report)`` or raises :class:`ConversionError`. +""" + +from __future__ import annotations + +import argparse +import json +import sys +from pathlib import Path +from typing import Any + +import gemmi +from pydantic import ValidationError + +from powderline.gateways.gsasii.schema import AXIAL_DIVERGENCE_MIN, gsasii_cell_groups, validate_recipe +from powderline.symmetry import SPECIAL_POSITION_TOL, UIJ_TOL, analyze_site, cell_tie_groups, coupling_groups + +SCHEMA_NAMES = {"GSASII_Rietveld": "gsasii.rietveld", "GSASII_SPF": "gsasii.spf"} +CELL = ("a", "b", "c", "alpha", "beta", "gamma") +UIJ = ("U11", "U22", "U33", "U12", "U13", "U23") +#: Iparm1 key -> (block, field) in the gsasii instrument (A103) and the 0.26.0 parameterization path. +IPARM1 = { + "Lam": ("radiation", "wavelength", ("wavelength",)), + "Polariz.": ("corrections", "polarization", ("polarization",)), + "Zero": ("corrections", "zero_shift", ("corrections", "zero_shift")), + "SH/L": ("corrections", "axial_divergence", ("corrections", "axial_divergence")), + **{k: ("broadening", k, ("broadening", k)) for k in "UVWXYZ"}, +} +IPARM1_PLAIN = ("Type", "Bank", "Azimuth") +PEAK_LISTS_026 = ("positions", "intensities", "pv_gaussian_sigma", "pv_lorentzian_gamma") + + +class Each: + """Spec of a name-keyed dict (phases, atoms): every value has ``spec``.""" + + def __init__(self, spec): + self.spec = spec + + +ANY = object() # carried whole (xrd_data extras are moved to metadata; Iparm1/Iparm2 are read by key) +UNUSED = object() # a 0.26.0 field the conversion does not carry (reported when it holds a value) +_FLAT = dict.fromkeys +_UIJ_SPEC = _FLAT(UIJ) +_SIZE_UNUSED = ("uniaxial_equatorial", "uniaxial_axial", "hkl_direction", "S11", "S22", "S33", "S12", "S13", "S23") +_STRAIN_UNUSED = ("uniaxial_equatorial", "uniaxial_axial", "hkl_direction", "stephens_parameters") +#: The 0.26.0 schema's fields (``powderline.schema``), for the unknown-key report (A132). +SPEC_026 = { + "schema_name": None, "schema_version": None, + "payload": { + "xrd_data": ANY, "asset_path": None, "fit_range": None, "phases": None, "single_peaks": None, + "instrument": {"description": None, "initialization": ANY, "parameterization": { + "wavelength": None, "polarization": None, "broadening": _FLAT("UVWXYZ"), + "corrections": _FLAT(("zero_shift", "axial_divergence"))}}, + "background": {"chebyshev": _FLAT(("num_coefficients", "coefficients", "refine_flag")), + "single_peaks": _FLAT(PEAK_LISTS_026)}, + "refinement_controls": {"refinement_cycles": None, "refinement_algorithm": None, + "single_peak_fitting_mode": _FLAT(("use_instrument_profile",))}, + }, +} +SPEC_026["payload"]["single_peaks"] = _FLAT(("positions", "intensities", "pv_gaussian_sigma_sq", "pv_lorentzian_gamma")) +SPEC_026["payload"]["phases"] = Each({ + "structure": {"phase_name": None, "space_group": None, "unit_cell": _FLAT(CELL + ("volume",)), + "atoms": Each({**_FLAT(("element", "x", "y", "z", "occupancy", "Multiplicity", "ADP", "Uiso")), + "Uaniso": _UIJ_SPEC, "Uij": _UIJ_SPEC})}, + "parameterization": { + "scale": None, "unit_cell": _FLAT(CELL), + "atoms": Each({**_FLAT(("x", "y", "z", "occupancy", "ADP", "Uiso")), "Uaniso": _UIJ_SPEC}), + "peak_broadening": { + "size_broadening": {"model": None, "isotropic_size": None, "LG_eta": None, + **dict.fromkeys(_SIZE_UNUSED, UNUSED)}, + "strain_broadening": {"model": None, "isotropic_strain": None, "LG_eta": None, + **dict.fromkeys(_STRAIN_UNUSED, UNUSED)}}}, +}) + + +class ConversionError(ValueError): + """The recipe cannot be converted as it stands; ``problems`` lists every reason.""" + + def __init__(self, problems: list[str]): + self.problems = problems + super().__init__("cannot convert the recipe:\n" + "\n".join(f" - {p}" for p in problems)) + + +def _value(param) -> Any: + return None if param is None else param[0] + + +def _flag(param) -> bool: + return bool(param[1]) if param is not None and param[1] is not None else False + + +def canonical_space_group(symbol: str) -> str: + """The canonical (gemmi ``xhm``) name of a 0.26.0 symbol, as GSAS-II read it (A105, A132). + + GSAS-II reads a trailing ``" R"`` or ``":R"`` as rhombohedral axes and a bare + (or ``":H"``) rhombohedral symbol as hexagonal axes, uses origin choice 2 for + a two-origin group, and reads a short monoclinic symbol as b-unique (gemmi's + choice too). Raises ``ValueError`` if gemmi cannot read the symbol, and for an + explicit origin choice 1: v0.1.1's GSAS-II read such coordinates as origin 2 + (EB-03), so neither the stated structure nor the one v0.1.1 refined is known. + """ + text = symbol.strip() + rhombohedral_axes = text.endswith((" R", ":R")) + base, _, setting = (text[:-2], "", "") if rhombohedral_axes else text.partition(":") + if setting.strip() == "1": + raise ValueError( + f"space group {symbol!r} states origin choice 1, but v0.1.1's GSAS-II read the coordinates in origin " + "choice 2 (EB-03), so the converter cannot tell which structure is meant; transform the coordinates to " + "origin choice 2 and state ':2'") + sg = gemmi.find_spacegroup_by_name(base.strip()) + if sg is None: + raise ValueError(f"space group {symbol!r} is not a symbol gemmi reads") + if sg.ext in ("H", "R"): + name = f"{sg.hm}:{'R' if rhombohedral_axes else 'H'}" + elif sg.ext in ("1", "2"): + name = f"{sg.hm}:2" + else: + name = sg.hm + return gemmi.find_spacegroup_by_name(name).xhm() + + +class _Conversion: + def __init__(self) -> None: + self.report: list[str] = [] + self.problems: list[str] = [] + + def note(self, path: str, text: str) -> None: + self.report.append(f"{path}: {text}") + + # --- top level ----------------------------------------------------------- + + def recipe(self, old: dict) -> dict: + name = old.get("schema_name") + if name not in SCHEMA_NAMES: + raise ConversionError([f"schema_name: {name!r} is not a 0.26.0 GSAS-II workflow ({', '.join(SCHEMA_NAMES)})"]) + new_name = SCHEMA_NAMES[name] + payload = old.get("payload") or {} + out: dict[str, Any] = {"schema_name": new_name, "core_schema_version": "1.0.0", + "engine_schema_version": "1.0.0"} + self.note("schema_name", f"{name} -> {new_name}; schema_version {old.get('schema_version')!r} -> " + "core_schema_version 1.0.0, engine_schema_version 1.0.0") + if payload.get("asset_path") is not None: + self.note("payload.asset_path", f"dropped ({payload['asset_path']!r}; not a gateway input, A50)") + new: dict[str, Any] = {} + metadata: dict[str, Any] = {} + new["xrd_data"] = self.xrd_data(payload.get("xrd_data") or {}, metadata) + new["instrument"] = self.instrument(payload.get("instrument") or {}) + if payload.get("fit_range") not in (None, [None, None]): + new["fit_range"] = list(payload["fit_range"]) + background = self.background(payload.get("background")) + if background: + new["background"] = background + if new_name == "gsasii.rietveld": + self.unused_peak_list(payload.get("single_peaks")) + new["phases"] = {k: self.phase(k, v) for k, v in (payload.get("phases") or {}).items()} + else: + if payload.get("phases"): + self.problems.append("payload.phases: a GSASII_SPF recipe has no phases") + new["single_peaks"] = self.peak_list(payload.get("single_peaks") or {}) + new["refinement_controls"] = self.controls(payload.get("refinement_controls") or {}, new_name) + out["payload"] = new + if metadata: + out["metadata"] = metadata + self.dropped_bounds(payload) + self.unknown_keys(old, SPEC_026, "") + return out + + def unknown_keys(self, old, spec, path: str) -> None: + """Report every key the 0.26.0 schema did not define (``extra='allow'``: v0.1.1 ignored it) and every + defined field the conversion does not carry that holds a value (A132).""" + if not isinstance(old, dict) or spec is ANY or spec is None: + return + if isinstance(spec, Each): + for k, v in old.items(): + self.unknown_keys(v, spec.spec, f"{path}.{k}") + return + for k, v in old.items(): + where = f"{path}.{k}".lstrip(".") + if k not in spec: + if path.endswith("parameterization.corrections") and v is not None: + self.problems.append(f"{where}: not a 0.26.0 correction; v0.1.1 raised for it") + else: + self.note(where, "dropped: not a 0.26.0 field (v0.1.1 ignored it)") + elif spec[k] is UNUSED: + if v is not None: + self.note(where, f"dropped ({v!r}): the isotropic model v0.1.1 ran did not use it") + else: + self.unknown_keys(v, spec[k], where) + + def xrd_data(self, old: dict, metadata: dict) -> dict: + kept = ("tth", "Itth", "Itth_weights", "filename") + for key in old: + if key not in kept: + metadata.setdefault("xrd_data", {})[key] = old[key] + self.note(f"payload.xrd_data.{key}", f"moved to metadata.xrd_data.{key}") + return {k: old[k] for k in kept if k in old} + + # --- instrument (A99, A103) ------------------------------------------------ + + def instrument(self, old: dict) -> dict: + init = old.get("initialization") or [{}] + iparm1 = init[0] if init else {} + param = old.get("parameterization") or {} + if len(init) > 1 and init[1]: + self.note("payload.instrument.initialization[1]", "Iparm2 dropped (v0.1.1 never used it)") + extra = sorted(set(iparm1) - set(IPARM1) - set(IPARM1_PLAIN)) + missing = sorted((set(IPARM1) | set(IPARM1_PLAIN)) - set(iparm1)) + if extra: + self.problems.append(f"payload.instrument.initialization[0]: GSAS-II keys {extra} are not supported " + "(e.g. a Kα doublet); gsasii 1.0.0 is single-wavelength PXC only") + if missing: + self.problems.append(f"payload.instrument.initialization[0]: missing GSAS-II keys {missing}") + if extra or missing: + return {} + if iparm1["Type"][1] != "PXC": + self.problems.append(f"payload.instrument.initialization[0].Type: {iparm1['Type'][1]!r}; gsasii 1.0.0 " + "supports 'PXC' only") + new: dict[str, Any] = {} + if old.get("description") is not None: + new["description"] = old["description"] + new["radiation"] = {"type": "PXC"} + new["geometry"] = {"bank": int(iparm1["Bank"][1]), "azimuth": iparm1["Azimuth"][1]} + for key, (block, field, path) in IPARM1.items(): + entry = param + for part in path: + entry = (entry or {}).get(part) + value = entry[0] if entry is not None and entry[0] is not None else iparm1[key][1] + flag = bool(entry[1]) if entry is not None and entry[1] is not None else bool(iparm1[key][2]) + new.setdefault(block, {})[field] = [value, flag] + shl = new["corrections"]["axial_divergence"] + if shl[0] < AXIAL_DIVERGENCE_MIN: + self.note("payload.instrument.corrections.axial_divergence", + f"SH/L {shl[0]} -> {AXIAL_DIVERGENCE_MIN} (GSAS-II computed with {AXIAL_DIVERGENCE_MIN} " + "anyway, EB-37; A91, deliberate change)") + shl[0] = AXIAL_DIVERGENCE_MIN + return new + + # --- background and peak lists --------------------------------------------- + + def background(self, old) -> dict: + new: dict[str, Any] = {} + cheb = (old or {}).get("chebyshev") + if cheb is None: + self.note("payload.background.chebyshev", "left out: no background (one fixed 0.0 term), where v0.1.1 " + "used GSAS-II's constant 1.0 (D1, A119; deliberate change)") + else: + new["chebyshev"] = {k: cheb[k] for k in ("num_coefficients", "coefficients", "refine_flag")} + peaks = (old or {}).get("single_peaks") + if peaks is not None: + lists = [peaks.get(k) or [] for k in PEAK_LISTS_026] + kept, dropped = [], 0 + for i in range(max(map(len, lists))): + entries = [lst[i] if i < len(lst) else None for lst in lists] + if any(e is None or e[0] is None or e[1] is None for e in entries): + dropped += 1 + else: + kept.append([[e[0], bool(e[1])] for e in entries]) + if dropped: + self.note("payload.background.single_peaks", f"{dropped} peak(s) with a null entry dropped " + "(v0.1.1 dropped them silently, D4)") + if kept: + new["single_peaks"] = {k: [peak[j] for peak in kept] for j, k in + enumerate(("positions", "intensities", "pv_gaussian_sigma_sq", + "pv_lorentzian_gamma"))} + self.note("payload.background.single_peaks", "pv_gaussian_sigma renamed pv_gaussian_sigma_sq " + "(the value always was σ², A89)") + return new + + def unused_peak_list(self, old) -> None: + if old and any(e is not None and any(v is not None for v in e) + for k in ("positions", "intensities", "pv_gaussian_sigma_sq", "pv_lorentzian_gamma") + for e in (old.get(k) or [])): + self.note("payload.single_peaks", "dropped: Peak List peaks take no part in a Rietveld refinement") + + def peak_list(self, old: dict) -> dict: + return {k: [[e[0], _flag(e)] for e in (old.get(k) or [])] + for k in ("positions", "intensities", "pv_gaussian_sigma_sq", "pv_lorentzian_gamma")} + + def controls(self, old: dict, schema_name: str) -> dict: + new: dict[str, Any] = {} + if "refinement_cycles" in old: + new["refinement_cycles"] = old["refinement_cycles"] + if old.get("refinement_algorithm") is not None: + self.note("payload.refinement_controls.refinement_algorithm", + f"dropped ({old['refinement_algorithm']!r}; v0.1.1 ignored it, D17)") + if schema_name == "gsasii.spf" and old.get("single_peak_fitting_mode") is not None: + new["single_peak_fitting_mode"] = old["single_peak_fitting_mode"] + return new + + def dropped_bounds(self, payload: dict) -> None: + """Report every populated min/max (A84: gsasii 1.0.0 has no bounds).""" + def walk(obj, path): + if isinstance(obj, dict): + for k, v in obj.items(): + if k not in ("xrd_data", "initialization"): + walk(v, f"{path}.{k}") + elif isinstance(obj, list): + if len(obj) == 4 and isinstance(obj[1], bool) and (obj[2] is not None or obj[3] is not None): + self.note(path, f"bounds [{obj[2]}, {obj[3]}] dropped (gsasii 1.0.0 exposes no bounds, A84)") + else: + for i, v in enumerate(obj): + walk(v, f"{path}[{i}]") + walk(payload, "payload") + + # --- phases (A93-A96, A100, A101, A105) ---------------------------------- + + def phase(self, name: str, old: dict) -> dict: + path = f"payload.phases.{name}" + structure = old.get("structure") or {} + param = old.get("parameterization") or {} + if structure.get("phase_name") != name: + self.problems.append(f"{path}.structure.phase_name: {structure.get('phase_name')!r} differs from the " + f"phase key {name!r} (the key is the name, A100)") + try: + sg = canonical_space_group(structure.get("space_group", "")) + except ValueError as exc: + self.problems.append(f"{path}.structure.space_group: {exc}") + return {} + if sg != structure.get("space_group"): + self.note(f"{path}.space_group", f"{structure.get('space_group')!r} -> {sg!r} (canonical name, as " + "GSAS-II read it; A105)") + new: dict[str, Any] = {"space_group": sg} + new["unit_cell"] = self.cell(path, sg, structure.get("unit_cell") or {}, param.get("unit_cell") or {}) + new["atoms"] = {label: self.atom(f"{path}.atoms.{label}", sg, atom, (param.get("atoms") or {}).get(label)) + for label, atom in (structure.get("atoms") or {}).items()} + for label in sorted(set(param.get("atoms") or {}) - set(structure.get("atoms") or {})): + self.note(f"{path}.parameterization.atoms.{label}", "dropped: the structure has no atom of that label " + "(v0.1.1 skipped it)") + scale = param.get("scale") + if scale is None or scale[0] is None: + self.note(f"{path}.scale", "no value: [1.0, flag] (GSAS-II's default, D9)") + new["scale"] = [1.0 if _value(scale) is None else scale[0], _flag(scale)] + broadening = self.broadening(path, param.get("peak_broadening") or {}) + if broadening: + new["peak_broadening"] = broadening + return new + + def merged(self, path: str, structure_value, param) -> Any: + if param is not None and param[0] is not None: + if structure_value is not None and param[0] != structure_value: + self.note(path, f"{structure_value} -> {param[0]} (the parameterization value won in v0.1.1)") + return param[0] + return structure_value + + def expand(self, path: str, params: dict, ties, extra_groups=()) -> None: + """Group-OR flags within tie groups (and GSAS-II's coarser groups); fixed parameters fixed (A95).""" + for members in [g.members for g in ties.groups] + list(extra_groups): + if any(params[m][1] for m in members): + for m in members: + if not params[m][1]: + params[m][1] = True + self.note(f"{path}.{m}", "refine flag false -> true (refined together with " + f"{', '.join(members)} in v0.1.1)") + for m in ties.fixed: + if params[m][1]: + params[m][1] = False + self.note(f"{path}.{m}", "refine flag true -> false (fixed by symmetry; v0.1.1 ignored it)") + + def cell(self, path: str, sg: str, structure: dict, param: dict) -> dict: + cell = {k: [self.merged(f"{path}.unit_cell.{k}", structure.get(k), param.get(k)), _flag(param.get(k))] + for k in CELL} + if structure.get("volume") is not None: + self.note(f"{path}.unit_cell.volume", "dropped (derived from a..gamma, A75)") + self.expand(f"{path}.unit_cell", cell, cell_tie_groups(sg), gsasii_cell_groups(sg)) + return cell + + def atom(self, path: str, sg: str, structure: dict, param) -> dict: + param = param or {} + new: dict[str, Any] = {"element": structure.get("element")} + for k in ("x", "y", "z"): + new[k] = [self.merged(f"{path}.{k}", structure.get(k), param.get(k)), _flag(param.get(k))] + occupancy = self.merged(f"{path}.occupancy", structure.get("occupancy"), param.get("occupancy")) + if occupancy is None: + occupancy = 1.0 + self.note(f"{path}.occupancy", "absent: 1.0 (v0.1.1's default)") + new["occupancy"] = [occupancy, _flag(param.get("occupancy"))] + xyz = tuple(new[k][0] for k in ("x", "y", "z")) + try: + site = analyze_site(sg, xyz) + ties = coupling_groups(sg, xyz) + except ValueError as exc: # e.g. the ambiguous band (A73): stop, never fix + self.problems.append(f"{path}: {exc}") + return new + self.expand(path, new, ties.xyz) + for i, k in enumerate(("x", "y", "z")): # a fixed constant v0.1.1 snapped (A121) + exact = site.exact[i] + if exact is not None and exact.denominator in (1, 2, 4, 8) and new[k][0] != float(exact) \ + and abs(new[k][0] - float(exact)) <= SPECIAL_POSITION_TOL: + self.note(f"{path}.{k}", f"{new[k][0]} -> {float(exact)} (the exact special-position value; " + "v0.1.1 snapped it, validation reads constants exactly, A121)") + new[k][0] = float(exact) + if structure.get("Multiplicity") is not None: + if structure["Multiplicity"] == site.multiplicity: + new["Multiplicity"] = structure["Multiplicity"] + else: + self.note(f"{path}.Multiplicity", f"stated {structure['Multiplicity']}, derived " + f"{site.multiplicity}: dropped (A101)") + self.adp(path, structure, param, new, ties) + return new + + def adp(self, path: str, structure: dict, param: dict, new: dict, ties) -> None: + """ADP type, values and flags as v0.1.1 used them (re/04 subplan §3).""" + p_adp, p_uiso, p_uaniso = param.get("ADP"), param.get("Uiso"), param.get("Uaniso") or {} + supplied_aniso = all(_value(p_uaniso.get(k)) is not None for k in UIJ) + given = [k for k in UIJ if _value(p_uaniso.get(k)) is not None] + if p_adp == "Uaniso" and 0 < len(given) < len(UIJ): + self.problems.append(f"{path}.Uaniso: the parameterization gives {', '.join(given)} only; v0.1.1 raised " + "unless all six are given (so there is no v0.1.1 refinement to reproduce)") + if p_adp == "Uiso" and _value(p_uiso) is not None: + adp = "Uiso" + elif p_adp == "Uaniso" and supplied_aniso: + adp = "Uaniso" + else: + adp = structure.get("ADP") + if p_adp is not None and p_adp != adp: + self.note(f"{path}.ADP", f"{adp} (the structure's): the parameterization's {p_adp} supplied no values, " + f"so GSAS-II kept {adp} and its {p_adp} flags applied to it") + new["ADP"] = adp + if adp == "Uiso": + value = self.merged(f"{path}.Uiso", structure.get("Uiso"), p_uiso) + if value is None: + value = 0.0 + self.note(f"{path}.Uiso", "absent: 0.0 (v0.1.1's default)") + flag = _flag(p_uiso) if p_adp != "Uaniso" else any(_flag(p_uaniso.get(k)) for k in UIJ) + new["Uiso"] = [value, flag] + return + old_u = structure.get("Uaniso") or structure.get("Uij") or {} + if not structure.get("Uaniso") and structure.get("Uij"): + self.note(f"{path}.Uij", "read as Uaniso (v0.1.1's reader accepted the old name)") + uaniso = {} + for k in UIJ: + value = self.merged(f"{path}.Uaniso.{k}", old_u.get(k), p_uaniso.get(k)) + if value is None: + value = 0.0 + self.note(f"{path}.Uaniso.{k}", "absent: 0.0") + flag = _flag(p_uiso) if p_adp == "Uiso" else _flag(p_uaniso.get(k)) + uaniso[k] = [value, flag] + self.expand(f"{path}.Uaniso", uaniso, ties.uij) + for k in ties.uij.fixed: # a fixed Uij v0.1.1 snapped (A121) + if uaniso[k][0] != 0.0 and abs(uaniso[k][0]) <= UIJ_TOL: + self.note(f"{path}.Uaniso.{k}", f"{uaniso[k][0]} -> 0.0 (fixed by symmetry, A121)") + uaniso[k][0] = 0.0 + new["Uaniso"] = uaniso + + def broadening(self, path: str, old: dict) -> dict: + new: dict[str, Any] = {} + for part, field, default_026, default_new in ( + ("size_broadening", "isotropic_size", "1 um", "10 um"), + ("strain_broadening", "isotropic_strain", "1000 microstrain", "0")): + block = old.get(part) + if block is None or _value(block.get(field)) is None: + self.note(f"{path}.peak_broadening.{part}", f"left out: the documented default {default_new} " + f"applies, where v0.1.1 used GSAS-II's {default_026} " + "(A87, A107; deliberate change)") + continue + if block.get("model") != "isotropic": + self.problems.append(f"{path}.peak_broadening.{part}.model: {block.get('model')!r}; gsasii 1.0.0 " + "supports 'isotropic' only") + continue + eta = block.get("LG_eta") + if _value(eta) is None: + self.note(f"{path}.peak_broadening.{part}.LG_eta", "no value: [1.0, flag] (GSAS-II's default)") + new[part] = {"model": "isotropic", field: [block[field][0], _flag(block[field])], + "LG_eta": [1.0 if _value(eta) is None else eta[0], _flag(eta)]} + return new + + +def convert(recipe: dict) -> tuple[dict, list[str]]: + """Convert a 0.26.0 recipe dict; return ``(gsasii recipe, report)``. + + Raises :class:`ConversionError` listing every problem, including every error + the gsasii 1.0.0 schema reports for the converted recipe. + """ + run = _Conversion() + new = run.recipe(recipe) + if not run.problems: + try: + validate_recipe(new) + except ValidationError as exc: + run.problems += [f"{'.'.join(map(str, e['loc']))}: {e['msg']}" for e in exc.errors()] + if run.problems: + raise ConversionError(run.problems) + return new, run.report + + +def main(argv=None) -> int: + parser = argparse.ArgumentParser(description="Convert a 0.26.0 GSAS-II recipe to gsasii 1.0.0.") + parser.add_argument("input", type=Path, help="0.26.0 recipe (JSON)") + parser.add_argument("-o", "--output", type=Path, help="write the converted recipe here (JSON)") + args = parser.parse_args(argv) + try: + new, report = convert(json.loads(args.input.read_text(encoding="utf-8"))) + except ConversionError as exc: + print(exc, file=sys.stderr) + return 1 + print("\n".join(report)) + if args.output: + with args.output.open("w", encoding="utf-8", newline="\n") as f: + json.dump(new, f, indent=2, ensure_ascii=False) + f.write("\n") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/update_code_hash.py b/scripts/update_code_hash.py index f537829..976f396 100644 --- a/scripts/update_code_hash.py +++ b/scripts/update_code_hash.py @@ -19,6 +19,7 @@ "executors.py", "extractors.py", "constraints.py", + "_legacy_executor.py", # verbatim 0.26.0 copy (re/04, A39); removed in re/07 ) diff --git a/src/powderline/_code_hash.json b/src/powderline/_code_hash.json index f83bfc6..8256b4b 100644 --- a/src/powderline/_code_hash.json +++ b/src/powderline/_code_hash.json @@ -1,11 +1,12 @@ { "gsasii_gateway_hashes": { + "_legacy_executor.py": "b541c8646e4a67ad809696793ea09c27", "constraints.py": "dc7a28029d0a071cf35dd352e76eaee9", - "executors.py": "a59bf71961245fe6a30842b092b7bb11", - "extractors.py": "7e57b8b0e48ec677d9a90f0441d8cee9", - "helpers.py": "d288784c535123f95d5d2ccc5083e806", - "project.py": "f074993af296685f632db2cb94504486", - "setters.py": "f7e1082b26194ff2ed86e9bde5e4eab2" + "executors.py": "0761a600e29e5f2d314feffc167509ab", + "extractors.py": "30cc981566bc069f2ba309cfc1f6242c", + "helpers.py": "5829fa4cb24cc41a69f27271061926f0", + "project.py": "90c618bf9b4433f6b4b876a4266a1781", + "setters.py": "536dcc49658e5006639e479b220c4703" }, "schema_version": "0.26.0" } diff --git a/src/powderline/gateways/gsasii/_legacy_executor.py b/src/powderline/gateways/gsasii/_legacy_executor.py new file mode 100644 index 0000000..1c3defe --- /dev/null +++ b/src/powderline/gateways/gsasii/_legacy_executor.py @@ -0,0 +1,2187 @@ +"""GSAS-II gateway — transient legacy executor for 0.26.0 ``GSASII_*`` recipes (re/04, A39). + +Verbatim copies, taken at ``2810c82`` (before the re/04 input-side edits), of +every function the re/04 frozen-edits ledger touches, and of every function on +the 0.26.0 path from ``run_refinement`` that calls one of them, so a 0.26.0 +recipe runs only unedited code: the before/after-refactor comparison stays +exact (maintainer, ledger sign-off 2026-10-08). The functions imported below are +unedited and shared. Generated (devkit ``tasks/re04-frozen-edits.md``); deleted +with the 0.26.0 schema in re/07. + +Do not edit: the bodies below are byte-identical to their originals. +""" +from __future__ import annotations +import sys +import uuid +import yaml +import json +import argparse +import random as ran +import numpy as np +import pandas as pd +from typing import Any, Dict, List, Literal, Optional, Sequence, Union +from pathlib import Path +from pydantic import ValidationError +from GSASII import GSASIIscriptable as G2 # type: ignore +from GSASII import GSASIIlattice as G2lat +from GSASII import GSASIImapvars as G2mv +from GSASII import GSASIIobj as G2obj +from GSASII import GSASIIspc as G2spc +from GSASII import GSASIIElem as G2elem +from GSASII import GSASIIpwd as G2pwd +from powderline.schema import RecipeModel, RefinementControls +from powderline.gateways.gsasii.constraints import atom_refinement_plan, cell_refinement_plan +from powderline._status import CHECK, CROSS, INFO, WARN, emoji +from dataclasses import dataclass + +from powderline.gateways.gsasii.helpers import ( # unedited, shared + DEFAULT_HIST_SCALE_REFINE_FLAG, + DEFAULT_HIST_SCALE_VAL, + DEFAULT_SPF_GAMMA_MIN, + DEFAULT_SPF_SIGMA_MIN, + OUTPUT_NAMING, + calculate_peak_widths, + require_phase, +) +from powderline.gateways.gsasii.project import ( # unedited, shared + add_powder_histogram_from_arrays, +) +from powderline.gateways.gsasii.extractors import ( # unedited, shared + build_atom_name_mapping, + build_phase_name_mapping, + calculate_cell_esds_from_A_matrix, + extract_refined_params_from_project, + get_descriptive_param_name, + parse_parameter_associations, +) +from powderline.gateways.gsasii.executors import ( # unedited, shared + execute_rietveld_refinement, +) +from powderline.gateways.gsasii.setters import ( # unedited, shared + set_chebyshev_background, + set_fit_range_hist, + set_hist_scale, + set_refinement_cycles, + set_single_peaks, +) + + +# --- verbatim from gateways/gsasii/helpers.py at 2810c82 --- + +def unpack_refinement_parameter(param: list | None, param_name: str = "parameter") -> tuple[Any, Any, Any, Any]: + """ + Unpack refinement parameter from [value, refine_flag, min, max] format. + + This utility reduces boilerplate for the common pattern of unpacking refinement + parameters throughout the codebase. Use only in straightforward cases - complex + conditional logic should inline the unpacking for clarity. + + Args: + param: RefinementParameter list [value, refine_flag, min, max] or None + param_name: Name of parameter for error messages (optional) + + Returns: + Tuple of (value, refine_flag, min_val, max_val) + If param is None, returns (None, None, None, None) + + Raises: + ValueError: If param is neither list nor None + + Examples: + >>> value, refine_flag, min_val, max_val = unpack_refinement_parameter([1.0, True, None, None]) + >>> value, refine_flag, min_val, max_val = unpack_refinement_parameter(None) + """ + if param is None: + return (None, None, None, None) + elif isinstance(param, list): + if len(param) != 4: + raise ValueError( + f"{param_name} must be a list of exactly 4 elements " + f"[value, refine_flag, min, max], got {len(param)} elements" + ) + return tuple(param) + else: + raise ValueError( + f"{param_name} must be a list [value, refine_flag, min, max] or None, " + f"got {type(param).__name__}" + ) + + +# --- verbatim from gateways/gsasii/helpers.py at 2810c82 --- + +def has_active_refinement_parameter(atom_param: dict) -> bool: + """ + Check if any refinement parameter has non-default values. + + Parameters: + atom_param : dict + Atom parameters dict with keys like 'x', 'y', 'z', 'occupancy', 'ADP', 'Uiso', 'Uaniso' + + Returns: + bool + True if any parameter should be set (not all [None, False, None, None]) + + Notes: + - Skips 'ADP' field (string, not refinement parameter) + - Handles 'Uaniso' nested dict separately + - Checks RefinementParameter lists for non-default values + """ + for key, value in atom_param.items(): + if key == 'ADP': + continue # Skip ADP type string (not a refinement parameter) + + if key == 'Uaniso' and isinstance(value, dict): + # Check nested Uaniso dict - any U value not [None, False, None, None]? + if any(v is not None and v != [None, False, None, None] for v in value.values()): + return True + + elif isinstance(value, list) and len(value) == 4: + # Check RefinementParameter list [value, refine_flag, min, max] + if value != [None, False, None, None]: + return True + + return False + + +# --- verbatim from gateways/gsasii/project.py at 2810c82 --- + +def add_phase_from_cif_dict(proj: Any, cif_data: Dict[str, Any], phase_name: str, histograms: Union[Sequence, str] = 'all') -> Any: + """ + Add a phase to a GSAS-II project from an embedded structure dictionary. + + This function creates a phase from an in-memory structure data dictionary + (file-less payload format) without reading from files. It creates a new phase entry + in the project, sets up the phase's general data (space group, unit_cell) and + atoms, and links the phase to specified histograms. + + Parameters + ---------- + proj : GSASIIscriptable.G2Project + The GSAS-II project object to which the phase will be added. + cif_data : dict + The structure dictionary containing phase structural information. + From the payload this dict is stowed in payload['phases'][phase]['structure'] + (e.g., payload['phases']['LaB6']['structure']). + phase_name : str + Name for the new phase. Top level key for a phase in payload['phases'], + or within payload['phases'][phase]['structure']['phase_name'] from payload. + histograms : list or 'all', optional + Which histograms to associate with this phase. Can be: + - 'all' (default): link to all powder histograms in the project. + - a list of histogram identifiers (names, objects, or indices). + - None or empty list for no links. + + Returns + ------- + phase_obj : GSASIIscriptable.G2Phase + The newly created phase object (GSAS-II scriptable phase wrapper). + + Raises + ------ + ValueError + If the space group in cif_data is invalid or cannot be interpreted. + + Notes + ----- + The function generates a new phase dictionary using GSASIIobj.SetNewPhase with + the provided space group and unit_cell parameters. Atomic coordinates and + occupancies from cif_data are inserted into the phase's atom list. Anisotropic + displacement parameters, if present, are fully supported. + + By default, the phase will be linked to all existing powder histograms in the + project. You can specify a subset of histograms by name, index, or object, or + pass an empty list/None to link none. + """ + pname = phase_name if phase_name else cif_data.get("name", "New Phase") + # Ensure phase name is unique in project + # TODO: handle phase name collisions more gracefully upstream in validation + existing_names = [p.name for p in proj.phases()] + pname = G2obj.MakeUniqueLabel(pname, existing_names) + + # Interpret space group symbol into SGData + sg = cif_data.get("space_group", "P 1") + err, SGData = G2spc.SpcGroup(sg) + + # Attempt normalization if initial interpretation failed + normalized_sg = sg + if err and sg: + normalized_sg = G2spc.StandardizeSpcName(sg) + if normalized_sg and normalized_sg != sg: + err, SGData = G2spc.SpcGroup(normalized_sg) + + # Handle space group interpretation errors + if err: + error_msg = G2spc.SGErrors(err) + raise ValueError( + f"Space group '{sg}' could not be interpreted.\n" + f"Error: {error_msg}\n" + f"Attempted normalization: '{normalized_sg}'\n" + f"\nCommon issues:\n" + f" • Missing spaces in Hermann-Mauguin notation (e.g., use 'R 3 m' not 'R3m')\n" + f" • Ambiguous rhombohedral symbols (R-3m normalizes to R-3c rhombohedral)\n" + f" • Check structure dict space_group entry matches the crystal structure\n" + ) + + # Create a new phase data structure (similar to reading from file): + phase_data = G2obj.SetNewPhase(Name=pname, SGData=SGData) + phase_data['General']['Name'] = pname + + # Document any space group normalization that occurred (silently - verbose mode not available here) + # if normalized_sg != sg: + # print(f"ℹ️ Space group normalization: '{sg}' → '{normalized_sg}'") + + # Set unit cell parameters if available + unit_cell = cif_data.get("unit_cell", {}) + if unit_cell: + unit_cell_constants = [unit_cell.get("a"), unit_cell.get("b"), unit_cell.get("c"), unit_cell.get("alpha"), unit_cell.get("beta"), unit_cell.get("gamma")] + if len(unit_cell_constants) == 6 and None not in unit_cell_constants: + phase_data['General']['Cell'][1:7] = unit_cell_constants # insert a,b,c,alpha,beta,gamma + try: + # Recalculate volume for consistency + phase_data['General']['Cell'][7] = G2lat.calc_V(G2lat.cell2A(tuple(unit_cell_constants))) + except Exception as e: + raise RuntimeError( + f"Unit cell volume calculation failed for phase '{pname}'. " + f"Cell parameters: {unit_cell_constants}. " + f"Original error: {e}" + ) from e + + # Add atoms from cif_data into phase_data['Atoms'] + # TODO: review handling of site_sym, multiplicity, and ADPS. + # The current approach expects minimal info in cif_data, but we are now populating this with the complete info needed. + # For example, cif_data['atoms'] stows info as shown below, where label is the key for each atom entry.: + # atoms[label] = { + # "element": element, + # "Multiplicity": mult, + # "x": x, "y": y, "z": z, + # "occupancy": occ, + # "ADP": ADP, + # "Uiso": Uiso_val + # } + # # Add Uij if anisotropic and data available + # if ADP == "Uani" and label in aniso_data: + # atoms[label]["Uij"] = aniso_data[label] + # TODO: in phase 2, assembly of atom_record using GSAS-II helpers will done upstream using "parse_cif_to_dict" function. + # This code will need to be update accordingly and use "get" with defaults for missing values. E.g., don't interpret site sym or mult. + + for label, atom in cif_data.get("atoms", {}).items(): + x = float(atom.get("x", 0.0)) + y = float(atom.get("y", 0.0)) + z = float(atom.get("z", 0.0)) + occ = float(atom.get("occupancy", 1.0)) + + # Get ADP type for this atom (required in schema 0.22) + adp_type = atom.get("ADP", None) + if adp_type is None: + raise ValueError( + f"Atom '{label}' in phase '{pname}' is missing required 'ADP' field. " + f"Schema 0.22 requires explicit ADP specification ('Uiso' or 'Uaniso') for all atoms." + ) + + # Prepare atom record list (based on GSAS-II internal format): + atom_record = ["", "", "", 0.0, 0.0, 0.0, 1.0, "", 0.0, "I", 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + atom_record[0] = label + atom_record[1] = atom.get("element", label.rstrip('0123456789') or "X") + atom_record[2] = "" # refinement flags placeholder (left blank) + atom_record[3] = x + atom_record[4] = y + atom_record[5] = z + atom_record[6] = occ + # Determine site symmetry and multiplicity for this atom in the given SG + site_sym, mult = "", 1 + try: + site_sym, mult = G2spc.SytSym([x, y, z], SGData)[0:2] + except Exception: + pass + atom_record[7] = site_sym if site_sym is not None else "" + atom_record[8] = mult if mult is not None else 1 + + # Set displacement parameters (isotropic or anisotropic) + if adp_type == "Uaniso": + # Anisotropic displacement parameters + atom_record[9] = "A" + + # Get Uaniso dict (supports both "Uij" old name and "Uaniso" new name) + uaniso_dict = atom.get("Uaniso") or atom.get("Uij") + + if uaniso_dict: + # Validate all 6 anisotropic values are present + required_keys = ["U11", "U22", "U33", "U12", "U13", "U23"] + missing_keys = [k for k in required_keys if k not in uaniso_dict or uaniso_dict[k] is None] + if missing_keys: + raise ValueError( + f"Atom '{label}' marked as anisotropic (ADP='Uaniso') but missing required U values: {missing_keys}. " + f"All 6 anisotropic parameters (U11, U22, U33, U12, U13, U23) must be provided." + ) + + # Populate U11, U22, U33, U12, U13, U23 (indices 11-16) + atom_record[11] = float(uaniso_dict["U11"]) + atom_record[12] = float(uaniso_dict["U22"]) + atom_record[13] = float(uaniso_dict["U33"]) + atom_record[14] = float(uaniso_dict["U12"]) + atom_record[15] = float(uaniso_dict["U13"]) + atom_record[16] = float(uaniso_dict["U23"]) + + # Set atom_record[10] to Uiso if provided, else calculate from diagonal elements + if atom.get("Uiso") is not None: + atom_record[10] = float(atom["Uiso"]) + else: + # Calculate equivalent isotropic value from diagonal elements + atom_record[10] = (atom_record[11] + atom_record[12] + atom_record[13]) / 3.0 + else: + raise ValueError( + f"Atom '{label}' marked as anisotropic (ADP='Uaniso') but no Uaniso dict found. " + f"Provide 'Uaniso' dict with U11, U22, U33, U12, U13, U23 values." + ) + else: + # Isotropic displacement parameters + atom_record[9] = "I" + uiso = atom.get("Uiso") + if uiso is not None: + atom_record[10] = float(uiso) + else: + atom_record[10] = 0.0 + atom_record[11:17] = [0.0]*6 # six anisotropic Uij (zeros since not used here) + + # Append a random unique ID for the atom + atom_record.append(ran.randint(0, sys.maxsize)) + phase_data['Atoms'].append(atom_record) + + # Insert the new phase into project data structure + if 'Phases' not in proj.data: + proj.data['Phases'] = {'data': None} + assert pname not in proj.data['Phases'], "Phase name collision despite uniqueness check" + proj.data['Phases'][pname] = phase_data + # Update the project tree name list for phases + for entry in proj.names: + if entry[0] == 'Phases': + entry.append(pname); break + else: + proj.names.append(['Phases', pname]) + + # Initialize phase general data (e.g. set default atom form factors): + try: + G2elem.SetupGeneral(phase_data, None) + except ValueError as err: + raise ValueError(f"Error in phase initialization: {err}") + + # Link phase to specified histograms + hist_list = [] + if histograms == 'all': + hist_list = [h.name for h in proj.histograms()] + elif histograms: + for h in histograms: + if hasattr(h, 'name'): + hist_list.append(h.name) + elif isinstance(h, str): + hist_list.append(h) + elif isinstance(h, int): + try: + hist_list.append(proj.histogram(h).name) + except Exception: + continue + for hist_name in hist_list: + try: + proj.link_histogram_phase(hist_name, pname) + except Exception as e: + print(f"Warning: could not link phase to histogram '{hist_name}': {e}") + + # Refresh internal IDs and return the new phase object + proj.index_ids() + proj.update_ids() + return proj.phase(pname) + + +# --- verbatim from gateways/gsasii/project.py at 2810c82 --- + +def add_phases_from_dict(proj: Any, hist: Any, phases_dict: dict, print_info: bool = False) -> None: + """ + Add phases to the GSAS-II project from a phases dictionary. + + Parameters: + proj : GSAS-II project object + The GSAS-II project to add phases to. + hist : GSAS-II histogram object + The histogram to associate with the phases. + phases_dict : dict + Dictionary containing phase information. (e.g., payload['phases']) + print_info : bool + Whether to print information about the added phases. + + Returns: + None + """ + + for phase_name, phase_info in phases_dict.items(): + structure_info = phase_info.get('structure', {}) + + # Validation check (TODO: move this upstream in phase2) + if phase_name != structure_info.get('phase_name'): + raise NameError(f"Phase name mismatch for '{phase_name}'. Structure info name is '{structure_info.get('phase_name')}'.") + + # Prevent duplicate phase names - TODO: handle upstream in schema validation / recipe maker + if phase_name in [p.name for p in proj.phases()]: + raise ValueError(f"Phase '{phase_name}' already exists in project.") + + try: + add_phase_from_cif_dict(proj=proj, cif_data=structure_info, phase_name=phase_name, histograms=[hist]) + except Exception as e: + raise RuntimeError(f"Failed to add phase '{phase_name}': {e}") from e + if print_info: + print("\nPhases in proj: ", *[ph.name for ph in proj.phases()], sep='\n\t') + + +# --- verbatim from gateways/gsasii/setters.py at 2810c82 --- + +def set_single_peak_background(proj: Any, hist: Any, bkg_single_peaks_dict: dict, print_info: bool = False) -> None: + """ + Set single peak background for histogram using pseudo-Voigt profiles. + + Single peaks are useful for modeling known impurity peaks or other non-background + features that shouldn't be included in the main phase refinement. Each peak is + described by a pseudo-Voigt profile (weighted sum of Gaussian and Lorentzian). + + Peak profile: I(2θ) = intensity * [η*L(2θ) + (1-η)*G(2θ)] + where: + - G(2θ) is Gaussian with width sigma + - L(2θ) is Lorentzian with width gamma + - η (eta) is mixing parameter (0=pure Gaussian, 1=pure Lorentzian) + + Args: + proj: GSAS-II project object + hist: Histogram object to set single peaks for + bkg_single_peaks_dict: Dictionary with keys: + - positions: List of [[2θ, refine, min, max], ...] for peak positions + - intensities: List of [[I, refine, min, max], ...] for peak heights + - pv_gaussian_sigma: List of [[σ, refine, min, max], ...] for Gaussian widths + - pv_lorentzian_gamma: List of [[γ, refine, min, max], ...] for Lorentzian widths + All lists must have same length (number of peaks) + print_info: If True, print peak configuration to stdout + + Returns: + None + + Raises: + ValueError: If parameter lists have inconsistent lengths + + Examples: + >>> # Two single peaks at 2θ=35.5° and 42.0° + >>> single_peaks = { + ... 'positions': [[35.5, False, None, None], [42.0, False, None, None]], + ... 'intensities': [[50.0, True, None, None], [30.0, True, None, None]], + ... 'pv_gaussian_sigma': [[0.1, False, None, None], [0.1, False, None, None]], + ... 'pv_lorentzian_gamma': [[0.05, False, None, None], [0.05, False, None, None]] + ... } + >>> set_single_peak_background(proj, hist, single_peaks) + """ + """ + Set the single peak background for a given histogram. + + Parameters: + proj: object + The GSAS-II project containing the histogram. + hist : object + The histogram to set the background for. + bkg_single_peaks_dict : dict + Dictionary containing single peak background parameters: + - positions : list of [float, bool, float, float]. Indices 2 and 3 can be None. + List of [position value, refine flag, min, max] for each peak. + These are two-theta values. Q / d can be passed upstream and converted in validation. + - intensities: list of [float, bool, float, float]. Indices 2 and 3 can be None. + List of [intensity value, refine flag, min, max] for each peak. + - pv_gaussian_sigma: list of [float, bool, float, float]. Indices 2 and 3 can be None. + List of [gaussian sigma value, refine flag, min, max] for each peak. + - pv_lorentzian_gamma: list of [float, bool, float, float]. Indices 2 and 3 can be None. + List of [lorentzian gamma value, refine flag, min, max] for each peak. + print_info : bool + Whether to print information about the set background. + Returns: + None + """ + # Parse the bkg_single_peaks_dict + positions = bkg_single_peaks_dict.get('positions', []) + intensities = bkg_single_peaks_dict.get('intensities', []) + gaussian_sigmas = bkg_single_peaks_dict.get('pv_gaussian_sigma', []) + lorentzian_gammas = bkg_single_peaks_dict.get('pv_lorentzian_gamma', []) + + # Default intensity, sigma, gamma values to use if value input is None + # TODO: Phase 2 upstream recipe builder — these defaults will be set before calling this function + # default_intensity = 1.0 + # default_sigma = 0.1 + # default_gamma = 0.1 + # default_refine_flag = False + + # Validate that all input lists have the same length + if not (len(positions) == len(intensities) == len(gaussian_sigmas) == len(lorentzian_gammas)): + raise ValueError("All input lists must have the same length.") + + num_peaks = len(positions) # only used for looping + + # Construct the single peak background list + single_peak_bkg = [] + + for i in range(num_peaks): + pos, pos_refine, pos_min, pos_max = positions[i] + inten, inten_refine, inten_min, inten_max = intensities[i] + g_sigma, g_refine, g_sigma_min, g_sigma_max = gaussian_sigmas[i] + l_gamma, l_refine, l_gamma_min, l_gamma_max = lorentzian_gammas[i] + + # TODO: move default behavior upstream + # # Use default values if None provided + # if inten is None: + # inten = default_intensity + # if g_sigma is None: + # g_sigma = default_sigma + # if l_gamma is None: + # l_gamma = default_gamma + + # # Use default refine flags if None provided + # if pos_refine is None: + # pos_refine = default_refine_flag + # if inten_refine is None: + # inten_refine = default_refine_flag + # if g_refine is None: + # g_refine = default_refine_flag + # if l_refine is None: + # l_refine = default_refine_flag + + # Min / max not currently available in GSAS-II single peak background definition + # They can be stored upstream for reference, but not used downstream currently. + # [placeholder for code to handle min/max if needed in future] + + # Append peak parameters to the single peak background list + # Each peak is a list: + # [pos, pos_refine, inten, inten_refine, g_sigma, g_refine, l_gamma, l_refine] + if None in [pos, pos_refine, + inten, inten_refine, + g_sigma, g_refine, + l_gamma, l_refine]: + continue + + # We only append if there are no Nones - this is a temp fix before phase2 + # TODO: update this behavior in the future when this list will only be passed with valid contents + single_peak_bkg.append([pos, pos_refine, + inten, inten_refine, + g_sigma, g_refine, + l_gamma, l_refine]) + + # Append the single peak background to the histogram's background list + if proj.data[hist.name]['Background'][1]['peaksList'] is None: + proj.data[hist.name]['Background'][1]['peaksList'] = [] + + # Since single_peak_bkg is a list of peaks, we need to append each peak individually + for peak in single_peak_bkg: + proj.data[hist.name]['Background'][1]['peaksList'].append(peak) + + # Update the number of peaks + proj.data[hist.name]['Background'][1]['nPeaks'] = len(proj.data[hist.name]['Background'][1]['peaksList']) + + # if print_info, then print details about the set background using the proj data, not the inputs + if print_info: + print(f"Single peak background set with {proj.data[hist.name]['Background'][1]['nPeaks']} peaks:") + for peak in proj.data[hist.name]['Background'][1]['peaksList']: + print(f" Position: {peak[0]} (refine: {peak[1]}), Intensity: {peak[2]} (refine: {peak[3]}), " + f"Gaussian Sigma: {peak[4]} (refine: {peak[5]}), Lorentzian Gamma: {peak[6]} (refine: {peak[7]})") + + +# --- verbatim from gateways/gsasii/setters.py at 2810c82 --- + +@require_phase("Cannot set scale.") +def set_phase_scale(proj: Any, hist: Any, phase_name: str, scale_param: list, print_info: bool = False) -> None: + """ + Set the phase scale factor in the GSAS-II project. + + Parameters: + proj : GSAS-II project object + The GSAS-II project containing the phase. + hist : GSAS-II histogram object + The histogram associated with the phase. + phase_name : str + The name of the phase to set the scale for. + scale_param : list + List containing [value, refine_flag, min, max] for the scale factor. + + Returns: + None + """ + value, refine_flag, min_val, max_val = unpack_refinement_parameter(scale_param, "scale_param") + + # Set the scale factor value - defaults will be handled upstream in future in validation + if value is not None: + proj.data['Phases'][phase_name]['Histograms'][hist.name]['Scale'][0] = value + + # Set refinement flag + if refine_flag is not None: + proj.data['Phases'][phase_name]['Histograms'][hist.name]['Scale'][1] = refine_flag + # Set min and max if provided - CURRENTLY NOT IMPLEMENETED. Unclear if allowed in GSAS-II. + # if min_val is not None: + # phase.set_parameter_min('Scale', min_val) + # if max_val is not None: + # phase.set_parameter_max('Scale', max_val) + + if print_info: + print(f"Set scale for phase '{phase_name}' to {value} with refine_flag={refine_flag}") + + +# --- verbatim from gateways/gsasii/setters.py at 2810c82 --- + +@require_phase("Cannot set unit cell parameters.") +def set_phase_unit_cell(proj: Any, phase_name: str, unit_cell_dict: dict, print_info: bool = False) -> list[str]: + """ + Set the unit cell parameters in the GSAS-II project for a given phase. + + Schema 0.26 semantics: each cell parameter refines iff it is present with + refine_flag=true; absent or false means fixed. Because GSAS-II has a single + whole-cell refine flag, per-parameter control is achieved by setting that + flag when any symmetry degree-of-freedom (DOF) group is requested and + returning "Hold" constraint variable names (e.g. '0::A2') for the DOF + groups that were not requested. Symmetry-linked parameters (e.g. cubic + a=b=c, or the coupled monoclinic {a, c, beta}) refine together if any + member is requested — see powderline.constraints. + + Parameters: + proj : GSAS-II project object + The GSAS-II project containing the phase. + phase_name : str + The name of the phase to set the unit cell for. + unit_cell_dict : dict + Dictionary containing unit cell parameters and their parameterization. + print_info : bool + Whether to print information about the set unit cell. + + Returns: + list[str] + GSAS-II variable names to Hold (empty when the whole cell refines or + is entirely fixed). The caller applies them via proj.add_HoldConstr. + """ + + # Note: cell values are set verbatim; validating them against the phase's + # crystal system is intentionally NOT done here (PowderLine is the engine, + # not the arbiter of recipe correctness — see docs/known_issues.md, KI-01). + + # Note2: hist not required since unit cell is phase-specific, not histogram-specific. + # The delineation of phase and histogram parameters is confusing in GSAS-II documentation. + # This makes sticking with a project based approach more straightforward. + + phase_names = [p.name for p in proj.phases()] + cell_changed = False + cell_params = ['a', 'b', 'c', 'alpha', 'beta', 'gamma'] # must match schema, and match order in GSAS-II + + for i, param in enumerate(cell_params): + if unit_cell_dict.get(param) is not None: + value, refine_flag, min_val, max_val = unit_cell_dict[param] + + # Set the unit cell parameter value - a=1, b=2, c=3, alpha=4, beta=5, gamma=6 + if value is not None: + proj.data['Phases'][phase_name]['General']['Cell'][i+1] = value + cell_changed = True + + # Set min and max if provided - CURRENTLY NOT IMPLEMENTED (deferred; + # GSAS-II's own min/max semantics differ from simple bounds). + + # Calculate and set volume from unit cell parameters only if cell was changed + if cell_changed: + cell = proj.data['Phases'][phase_name]['General']['Cell'][1:7] # get list of a, b, c, alpha, beta, gamma + proj.data['Phases'][phase_name]['General']['Cell'][7] = G2lat.calc_V(G2lat.cell2A(cell)) # index 7 is volume + + # Translate per-parameter refine flags into the whole-cell flag + holds + SGData = proj.data['Phases'][phase_name]['General']['SGData'] + phase_idx = phase_names.index(phase_name) + cell_plan = cell_refinement_plan(SGData, unit_cell_dict, phase_idx) + proj.data['Phases'][phase_name]['General']['Cell'][0] = cell_plan.refine_cell + + if print_info: + curr_cell = proj.data['Phases'][phase_name]['General']['Cell'] + if cell_changed: + cell_prms_dict = {"a": curr_cell[1], + "b": curr_cell[2], + "c": curr_cell[3], + "alpha": curr_cell[4], + "beta": curr_cell[5], + "gamma": curr_cell[6], + "volume": curr_cell[7]} + + print(f"Set unit cell parameters for phase '{phase_name}' to {[(key, value) for key, value in cell_prms_dict.items()]} with refine_flag={curr_cell[0]}") + else: + print(f"No unit cell parameters were changed for phase '{phase_name}'. Refine_flag={curr_cell[0]}") + if cell_plan.holds: + print(f"Holding fixed unit cell DOFs for phase '{phase_name}': {cell_plan.holds}") + + return cell_plan.holds + + +# --- verbatim from gateways/gsasii/setters.py at 2810c82 --- + +@require_phase("Cannot set size broadening parameters.") +def set_phase_size_broadening(proj: Any, hist: Any, phase_name: str, size_broadening_dict: dict, print_info: bool = False) -> None: + """ + Set size broadening parameters in the GSAS-II project for a given phase. + + Supports model branching for isotropic, uniaxial, and ellipsoidal models. + Currently only isotropic model is implemented - others raise NotImplementedError. + + Parameters: + proj : GSAS-II project object + The GSAS-II project containing the phase. + hist : GSAS-II histogram object + The histogram associated with the phase. + phase_name : str + The name of the phase to set the size broadening for. + size_broadening_dict : dict + Dictionary containing size broadening parameters with 'model' key. + + Raises: + NotImplementedError: For uniaxial or ellipsoidal models + ValueError: For unknown model types or invalid parameters + + Returns: + None + """ + # Extract model type (default to isotropic for backward compatibility) + model = size_broadening_dict.get('model', 'isotropic') + + # Validate model type + valid_models = ['isotropic', 'uniaxial', 'ellipsoidal'] + if model not in valid_models: + raise ValueError( + f"Unknown size broadening model '{model}'. " + f"Valid options: {valid_models}" + ) + + # Set model type in GSAS-II data structure + proj.data['Phases'][phase_name]['Histograms'][hist.name]['Size'][0] = model + + # Branch on model type + if model == 'isotropic': + _set_isotropic_size_broadening(proj, hist, phase_name, size_broadening_dict, print_info) + elif model == 'uniaxial': + raise NotImplementedError( + "Uniaxial size broadening model is not yet implemented. " + "Support planned for future release. Use 'isotropic' model instead." + ) + elif model == 'ellipsoidal': + raise NotImplementedError( + "Ellipsoidal size broadening model is not yet implemented. " + "Support planned for future release. Use 'isotropic' model instead." + ) + + +# --- verbatim from gateways/gsasii/setters.py at 2810c82 --- + +def _set_isotropic_size_broadening(proj: Any, hist: Any, phase_name: str, size_broadening_dict: dict, print_info: bool = False) -> None: + """ + Set isotropic size broadening parameters (internal helper). + + Refactored from original set_phase_size_broadening function. + """ + + # Example structure of proj.data['Phases'][phase_name]['Histograms'][hist.name]['Size']: + # ['isotropic', # [0]: str for broadening type. 'isotropic', 'uniaxial', or 'ellipsoidal' + # [1.0, 1.0, 1.0], # [1][0]: iso size or equatorial size, [1][1]: axial size (if uniaxial used, else ignored), [1][2]: LG_mix (eta parameter, 1 = Lorentzian, 0 = Gaussian) + # [False, False, False], # boolean refinement flags. [2][0]: isotropic or uniaxial equatorial refine, [2][1]: uniaxial axial refine, [2][2]: LG_mix refine + # [0, 0, 1], # hkl direction for uniaxial broadening. [3][0]: h, [3][1]: k, [3][2]: l + # [1.0, 1.0, 1.0, 0.0, 0.0, 0.0], # ellipsoidal sizes. S11, S22, S33, S12, S13, S23 for [4][0:6] + # [False, False, False, False, False, False]] # ellipsoidal size refine flags. S11, S22, S33, S12, S13, S23 for [5][0:6] + ############################ + # Set isotropic size parameters + ############################ + value, refine_flag, min_val, max_val = unpack_refinement_parameter( + size_broadening_dict.get('isotropic_size'), "isotropic_size" + ) + + # Set the size parameter value + if value is not None: + proj.data['Phases'][phase_name]['Histograms'][hist.name]['Size'][1][0] = value + + # Set refinement flag + if refine_flag is not None: + proj.data['Phases'][phase_name]['Histograms'][hist.name]['Size'][2][0] = refine_flag + + # Set min and max if provided - CURRENTLY NOT IMPLEMENTED. Unclear if allowed in GSAS-II. + + if print_info: + print(f"Set size broadening size parameters for phase '{phase_name}' to {value} with refine_flag={refine_flag}") + + ############################ + # Set size LG_eta parameter + ############################ + try: + value, refine_flag, min_val, max_val = unpack_refinement_parameter( + size_broadening_dict.get('LG_eta'), "LG_eta" + ) + except (ValueError, TypeError) as e: + # Preserve original error message format + raise ValueError(f"Size broadening 'LG_eta' parameter must be a list or None.\nType for phase {phase_name} is {type(size_broadening_dict['LG_eta'])}.") from e + + # Set the LG_eta parameter value + if value is not None: + proj.data['Phases'][phase_name]['Histograms'][hist.name]['Size'][1][2] = value + + # Set refinement flag + if refine_flag is not None: + proj.data['Phases'][phase_name]['Histograms'][hist.name]['Size'][2][2] = refine_flag + + # Set min and max if provided - CURRENTLY NOT IMPLEMENTED. Unclear if allowed in GSAS-II. + + if print_info: + print(f"Set size broadening LG_eta parameters for phase '{phase_name}' to {value} with refine_flag={refine_flag}") + + +# --- verbatim from gateways/gsasii/setters.py at 2810c82 --- + +@require_phase("Cannot set strain broadening parameters.") +def set_phase_strain_broadening(proj: Any, hist: Any, phase_name: str, strain_broadening_dict: dict, print_info: bool = False) -> None: + """ + Set strain broadening parameters in the GSAS-II project for a given phase. + + Supports model branching for isotropic, uniaxial, and generalized (Stephens) models. + Currently only isotropic model is implemented - others raise NotImplementedError. + + Parameters: + proj : GSAS-II project object + The GSAS-II project containing the phase. + hist : GSAS-II histogram object + The histogram associated with the phase. + phase_name : str + The name of the phase to set the strain broadening for. + strain_broadening_dict : dict + Dictionary containing strain broadening parameters with 'model' key. + + Raises: + NotImplementedError: For uniaxial or generalized models + ValueError: For unknown model types or invalid parameters + + Returns: + None + """ + # Extract model type (default to isotropic for backward compatibility) + model = strain_broadening_dict.get('model', 'isotropic') + + # Validate model type + valid_models = ['isotropic', 'uniaxial', 'generalized'] + if model not in valid_models: + raise ValueError( + f"Unknown strain broadening model '{model}'. " + f"Valid options: {valid_models}" + ) + + # Set model type in GSAS-II data structure + proj.data['Phases'][phase_name]['Histograms'][hist.name]['Mustrain'][0] = model + + # Branch on model type + if model == 'isotropic': + _set_isotropic_strain_broadening(proj, hist, phase_name, strain_broadening_dict, print_info) + elif model == 'uniaxial': + raise NotImplementedError( + "Uniaxial strain broadening model is not yet implemented. " + "Support planned for future release. Use 'isotropic' model instead." + ) + elif model == 'generalized': + raise NotImplementedError( + "Generalized (Stephens) strain broadening model is not yet implemented. " + "This model requires complex symmetry-dependent parameterization. " + "Support planned for Phase 2. Use 'isotropic' model instead." + ) + + +# --- verbatim from gateways/gsasii/setters.py at 2810c82 --- + +def _set_isotropic_strain_broadening(proj: Any, hist: Any, phase_name: str, strain_broadening_dict: dict, print_info: bool = False) -> None: + """ + Set isotropic strain broadening parameters (internal helper). + + Refactored from original set_phase_strain_broadening function. + """ + ############################ + # Set isotropic strain parameters + ############################ + value, refine_flag, min_val, max_val = unpack_refinement_parameter( + strain_broadening_dict.get('isotropic_strain'), "isotropic_strain" + ) + + # Set the strain parameter value + if value is not None: + proj.data['Phases'][phase_name]['Histograms'][hist.name]['Mustrain'][1][0] = value + + # Set refinement flag + if refine_flag is not None: + proj.data['Phases'][phase_name]['Histograms'][hist.name]['Mustrain'][2][0] = refine_flag + + # Set min and max if provided - CURRENTLY NOT IMPLEMENTED. Unclear if allowed in GSAS-II. + + if print_info: + print(f"Set strain broadening strain parameters for phase '{phase_name}' to {value} with refine_flag={refine_flag}") + + ############################ + # Set strain LG_eta parameter + ############################ + try: + value, refine_flag, min_val, max_val = unpack_refinement_parameter( + strain_broadening_dict.get('LG_eta'), "LG_eta" + ) + except (ValueError, TypeError) as e: + # Preserve original error message format + raise ValueError(f"Strain broadening dict key 'LG_eta' must hold a value that is a list or None.\nType for phase {phase_name} is {type(strain_broadening_dict['LG_eta'])}.") from e + + # Set the LG_eta parameter value + if value is not None: + proj.data['Phases'][phase_name]['Histograms'][hist.name]['Mustrain'][1][2] = value + + # Set refinement flag + if refine_flag is not None: + proj.data['Phases'][phase_name]['Histograms'][hist.name]['Mustrain'][2][2] = refine_flag + + # Set min and max if provided - CURRENTLY NOT IMPLEMENTED. Unclear if allowed in GSAS-II. + + if print_info: + print(f"Set strain broadening LG_eta parameters for phase '{phase_name}' to {value} with refine_flag={refine_flag}") + + +# --- verbatim from gateways/gsasii/setters.py at 2810c82 --- + +@require_phase("Cannot set atom parameters.") +def set_phase_atom_parameters(proj: Any, phase_name: str, atom_parameters_dict: dict, print_info: bool = False) -> list[str]: + """ + Set atom parameters for a phase in the GSAS-II project. + + This function parses atom-specific refinement parameters (coordinates, occupancy, displacement) + and constructs GSAS-II's refine_flags string ('F', 'X', 'U' combinations). + + Schema 0.26 semantics: each coordinate (x/y/z) and each anisotropic Uij + component refines iff present with refine_flag=true; absent or false means + fixed. GSAS-II's 'X'/'U' flags are per-atom, so per-component control is + achieved by returning "Hold" constraint variable names (e.g. '0::dAy:3', + '0::AU22:4') for the site-symmetry DOF groups that were not requested — + see powderline.constraints. + Symmetry-linked components (e.g. x=y on an (x,x,z) site) refine together + if any member is requested. + + Parameters: + proj : GSAS-II project object + The GSAS-II project containing the phase. + phase_name : str + The name of the phase to set the atom parameters for. + atom_parameters_dict : dict + Dictionary containing atom parameters and their parameterization. + E.g., payload['phases'][phase_name]['parameterization']['atoms'] + + Expected structure: + { + 'atom_label': { + 'x': [value, refine_flag, min, max], + 'y': [value, refine_flag, min, max], + 'z': [value, refine_flag, min, max], + 'occupancy': [value, refine_flag, min, max], + 'ADP': 'Uiso' or 'Uaniso', + 'Uiso': [value, refine_flag, min, max], # if ADP='Uiso' + 'Uaniso': { # if ADP='Uaniso' + 'U11': [value, refine_flag, min, max], + 'U22': [value, refine_flag, min, max], + ... + } + } + } + print_info : bool, optional + If True, print detailed information about the atom parameters being set. Default is False. + + Returns: + list[str] + GSAS-II variable names to Hold (empty when every requested DOF group + refines in full). The caller applies them via proj.add_HoldConstr. + + Notes: + proj.data['Phases'][phase_name]['Atoms'] is a list of lists, where each sublist represents an atom and its parameters. + Sublist structure in atom_param_index_mapping below. + + Example: proj.data['Phases'][phase_name]['Atoms'][0] from LaB6 example returns: + ['La', 'La', '', 0.0, 0.0, 0.0, 1.0, 'm3m', 1, 'I', 0.00858, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 844046454091148456] + + GSAS-II refine_flags string format (index 2): + - 'F': refine occupancy + - 'X': refine all coordinates (x, y, z together) + - 'U': refine displacement parameters (isotropic or anisotropic) + - Combinations: 'FXU', 'XU', 'F', etc. + """ + + # Dictionary of atom parameter names to their indices in the atom list + atom_param_index_mapping = { + 'label': 0, # atom label + 'element': 1, # element symbol + 'refine_flags': 2, # string of refine flags e.g., 'F', 'X', 'U' or any combination such as 'FXU'. F for occupancy, X for xyz coords, U for displacement parameters + 'x': 3, # fractional coordinate x value + 'y': 4, # fractional coordinate y value + 'z': 5, # fractional coordinate z value + 'occupancy': 6, # float, referred to as frac in GSAS-II syntax + 'site_symmetry': 7, # string representation of site symmetry + 'multiplicity': 8, # multiplicity of the site, integer + 'Utype': 9, # 'I' for isotropic, 'A' for anisotropic + 'Uiso': 10, # value if isotropic + 'U11': 11, # value if anisotropic + 'U22': 12, # value if anisotropic + 'U33': 13, # value if anisotropic + 'U12': 14, # value if anisotropic + 'U13': 15, # value if anisotropic + 'U23': 16 # value if anisotropic + # Note: index 17 is a random unique ID for the atom, not a parameter to set + } + + phase_names = [p.name for p in proj.phases()] + phase_idx = phase_names.index(phase_name) + SGData = proj.data['Phases'][phase_name]['General']['SGData'] + atom_labels_present = [x[0] for x in proj.data['Phases'][phase_name]['Atoms']] + holds: list[str] = [] + + for atom_label, atom_params in atom_parameters_dict.items(): + if atom_label not in atom_labels_present: + print(f"Atom '{atom_label}' not found in phase '{phase_name}'. Skipping.") + continue + + # Get the atom record + atom_index = atom_labels_present.index(atom_label) + atom_record = proj.data['Phases'][phase_name]['Atoms'][atom_index] + + # Get ADP type for this atom (required in schema 0.22) + adp_type = atom_params.get('ADP') + + # Process coordinate parameters (x, y, z) — values only; refine flags + # are translated by atom_refinement_plan below + for coord in ['x', 'y', 'z']: + if coord in atom_params and atom_params[coord] is not None: + param_value = atom_params[coord] + if param_value != [None, False, None, None]: + value, refine_flag, min_val, max_val = param_value + + # Update coordinate value if provided + if value is not None: + atom_record[atom_param_index_mapping[coord]] = value + + # Process occupancy + if 'occupancy' in atom_params and atom_params['occupancy'] is not None: + param_value = atom_params['occupancy'] + if param_value != [None, False, None, None]: + value, refine_flag, min_val, max_val = param_value + + # Update occupancy value if provided + if value is not None: + atom_record[6] = value + + # Process displacement parameters based on ADP type + if adp_type == 'Uiso': + # Isotropic displacement parameters + if 'Uiso' in atom_params and atom_params['Uiso'] is not None: + param_value = atom_params['Uiso'] + if param_value != [None, False, None, None]: + value, refine_flag, min_val, max_val = param_value + + # Update Uiso value if provided + if value is not None: + atom_record[9] = 'I' # Set to isotropic + atom_record[10] = value + + elif adp_type == 'Uaniso': + # Anisotropic displacement parameters + if 'Uaniso' in atom_params and atom_params['Uaniso'] is not None: + uaniso_dict = atom_params['Uaniso'] + + # Track if any anisotropic values are being updated + aniso_values_provided = [] + + # Process each anisotropic parameter + for u_key in ['U11', 'U22', 'U33', 'U12', 'U13', 'U23']: + if u_key in uaniso_dict and uaniso_dict[u_key] is not None: + param_value = uaniso_dict[u_key] + if param_value != [None, False, None, None]: + value, refine_flag, min_val, max_val = param_value + + # Update anisotropic value if provided + if value is not None: + atom_record[9] = 'A' # Set to anisotropic + atom_record[atom_param_index_mapping[u_key]] = value + aniso_values_provided.append(u_key) + + # Validate that if any anisotropic values are provided, all 6 must be provided + if aniso_values_provided: + required_keys = ['U11', 'U22', 'U33', 'U12', 'U13', 'U23'] + missing_keys = [k for k in required_keys if k not in aniso_values_provided] + if missing_keys: + raise ValueError( + f"Atom '{atom_label}' in phase '{phase_name}': " + f"When updating anisotropic displacement parameters, all 6 values (U11-U23) must be provided. " + f"Missing: {missing_keys}" + ) + + # Translate per-parameter refine flags into GSAS-II's per-atom flag + # string plus holds for the unrefined site-symmetry DOF groups. + # xyz/SGData let the plan recompute a stale/blank site symmetry the + # same way GSAS-II itself does (GetPhaseData's KeyError patch). + atom_plan = atom_refinement_plan( + atom_params, + atom_record[7], + phase_idx, + atom_index, + xyz=atom_record[3:6], + SGData=SGData, + ) + atom_record[2] = atom_plan.refine_flags + holds.extend(atom_plan.holds) + + # Verbose logging + if print_info and atom_plan.refine_flags: + adp_display = adp_type if adp_type else "inherited" + print(f" Set refine flags for atom '{atom_label}' ({adp_display}): '{atom_plan.refine_flags}'") + if atom_plan.holds: + print(f" Holding fixed DOFs for atom '{atom_label}': {atom_plan.holds}") + + return holds + + +# --- verbatim from gateways/gsasii/setters.py at 2810c82 --- + +def set_phase_parameterization(proj: Any, hist: Any, phases_dict: dict, print_info: bool = False) -> list[str]: + """ + Set phase parameterization in the GSAS-II project from a phases dictionary. + + Parameters: + proj : GSAS-II project object + The GSAS-II project containing the phases. + hist : GSAS-II histogram object + The histogram associated with the phases. + phases_dict : dict + Dictionary containing phase parameterization information. E.g., payload['phases'] + print_info : bool, optional + If True, print information about the parameterization being set. Default is False. + + Returns: + list[str] + Accumulated GSAS-II "Hold" variable names from the unit-cell and atom + setters (schema 0.26 per-parameter refine flags). The caller applies + them once via proj.add_HoldConstr. + """ + holds: list[str] = [] + + for phase_name, phase_info in phases_dict.items(): + + # In schema 0.21, phase_name is the dict key itself + # (structure.phase_name should match, but we use the key for consistency) + + # Check that phase exists in project - this would indicate an error in phase addition earlier + if phase_name not in [p.name for p in proj.phases()]: + print(f"Phase '{phase_name}' not found in project. Skipping parameterization.") + continue + + # Now that we have the phase name, we can add phase-specific parameterization here + # This will call on a few helper functions to set scale, unit cell, peak broadening, atom parameters, etc. + # Each of these helper functions will take the proj, phase name (rather than phase object), and relevant parameterization dict + + # First, get parameterization dict from phase_info (already phase-specific) + param_dict = phase_info.get('parameterization', {}) + + # Keys in param dict are 'scale', 'unit_cell', 'peak_broadening', and 'atoms' for now + # It is not as simple as looping over keys since some have nested dicts and others hold values or lists + + # Set scale factor if provided + scale_param = param_dict.get('scale', None) + if scale_param is not None: + if any(value is not None for value in scale_param): + set_phase_scale(proj, hist, phase_name, scale_param, print_info=False) + + # Set unit cell parameters if provided. The key values can be None, + # so here we cannot just check for unit_cell_dict not being None. + # We need to check if any of hte keys in unit_cell_dict are not None. + unit_cell_dict = param_dict.get('unit_cell', None) + + if unit_cell_dict is not None: + if any(value is not None for value in unit_cell_dict.values()): + holds.extend(set_phase_unit_cell(proj, phase_name, unit_cell_dict, print_info=False)) + + # Set peak broadening if provided. Keys are "size_broadening" and "strain_broadening" + # and each holds a dict with relevant parameters (size or strain) and the LG_eta + # ("or {}" because model_dump emits an explicit None when the section is absent) + peak_broadening_dict = param_dict.get('peak_broadening') or {} + size_broadening_dict = peak_broadening_dict.get('size_broadening', None) + strain_broadening_dict = peak_broadening_dict.get('strain_broadening', None) + + if size_broadening_dict is not None: + if any(value is not None for value in size_broadening_dict.values()): + set_phase_size_broadening(proj, hist, phase_name, size_broadening_dict, print_info=False) + + if strain_broadening_dict is not None: + if any(value is not None for value in strain_broadening_dict.values()): + set_phase_strain_broadening(proj, hist, phase_name, strain_broadening_dict, print_info=False) + + atoms_parameters_dict = param_dict.get('atoms', {}) + if atoms_parameters_dict and any(has_active_refinement_parameter(ap) for ap in atoms_parameters_dict.values()): + holds.extend(set_phase_atom_parameters(proj, phase_name, atoms_parameters_dict, print_info=False)) + + if print_info: + print(f"Completed parameterization for phase '{phase_name}'.") + + if print_info: + print("All phase parameterization complete. Phases in proj: ", *[p.name for p in proj.phases()], sep='\n\t') + if holds: + print("Hold constraints to apply (per-parameter refine flags): ", *holds, sep='\n\t') + + return holds + + +# --- verbatim from gateways/gsasii/setters.py at 2810c82 --- + +def set_instrument_parameterization(proj: Any, hist: Any, instrument_param_dict: dict, print_info: bool = False) -> None: + """ + Set instrument parameterization in the GSAS-II project from an instrument parameterization dictionary. + + Parameters: + proj : GSAS-II project object + The GSAS-II project containing the histogram. + hist : GSAS-II histogram object + The histogram associated with the instrument. + instrument_param_dict : dict + Dictionary containing instrument parameterization information. E.g., payload['instrument'] + print_info : bool, optional + If True, print information about the parameterization being set. Default is False. + + Returns: + None + """ + + # Note: GSAS-II's instrument parameters hold values in a list. + # The indices for each parameter (e.g., 'Lam' for wavelength) are as follows: + # [0] default/starting value from iprms, [1] = value, [2] = refinement flag + + # Wavelength + wavelength_param = instrument_param_dict.get('wavelength', None) + if wavelength_param is not None: + value, refine_flag, min_val, max_val = unpack_refinement_parameter(wavelength_param, "wavelength") + wavelength_changed = False # track if wavelength value changed + + # Set wavelength value + curr_wavelength = proj.data[hist.name]['Instrument Parameters'][0]['Lam'][1] + + if value is not None and value != curr_wavelength: + proj.data[hist.name]['Instrument Parameters'][0]['Lam'][1] = value + wavelength_changed = True + + # Set refinement flag + if refine_flag is not None: + proj.data[hist.name]['Instrument Parameters'][0]['Lam'][2] = refine_flag + + if print_info: + if wavelength_changed: + print(f"Set instrument wavelength to {value} with refine_flag={refine_flag}") + else: + print(f"No change to instrument wavelength. Refine_flag={refine_flag}") + + # Additional instrument parameters (polarization, broadening, corrections) would be set similarly + # Implementing those follows the same pattern as above + + # Polarization + polarization_param = instrument_param_dict.get('polarization', None) + if polarization_param is not None: + value, refine_flag, min_val, max_val = unpack_refinement_parameter(polarization_param, "polarization") + polarization_changed = False # track if polarization value changed + + curr_polarization = proj.data[hist.name]['Instrument Parameters'][0]['Polariz.'][1] + + # Set polarization value + if value is not None and value != curr_polarization: + proj.data[hist.name]['Instrument Parameters'][0]['Polariz.'][1] = value + polarization_changed = True + + # Set refinement flag + if refine_flag is not None: + proj.data[hist.name]['Instrument Parameters'][0]['Polariz.'][2] = refine_flag + + if print_info: + if polarization_changed: + print(f"Set instrument polarization to {value} with refine_flag={refine_flag}") + else: + print(f"No change to instrument polarization. Refine_flag={refine_flag}") + + # Broadening parameters (U, V, W, X, Y, Z) + broadening_dict = instrument_param_dict.get('broadening', {}) + for param_key in broadening_dict.keys(): + broadening_param = broadening_dict.get(param_key, None) + if broadening_param is not None: + value, refine_flag, min_val, max_val = unpack_refinement_parameter(broadening_param, f"broadening.{param_key}") + + broadening_changed = False + + # Set broadening parameter value + if value is not None: + proj.data[hist.name]['Instrument Parameters'][0][param_key][1] = value + broadening_changed = True + + # Set refinement flag + if refine_flag is not None: + proj.data[hist.name]['Instrument Parameters'][0][param_key][2] = refine_flag + + if print_info: + if broadening_changed: + # Print changed parameter + print(f"Set instrument broadening parameter '{param_key}' to {value} with refine_flag={refine_flag}") + + # Always print refinement flag even if value not changed + else: + print(f"No change to instrument broadening parameter '{param_key}'. Refine_flag={refine_flag}") + + + # Corrections parameters (zero_shift, axial_divergence, sample_height_displacement) + corrections_dict = instrument_param_dict.get('corrections', {}) + for param_key in corrections_dict.keys(): + correction_param = corrections_dict.get(param_key, None) + if correction_param is not None: # only loop over provided params + + if param_key == 'zero_shift': + gsas_key = 'Zero' + elif param_key == 'axial_divergence': + gsas_key = 'SH/L' + #elif param_key == 'sample_height_displacement': # irrelevant to area detector data, future work for Bragg-Brentano geometry + # gsas_key = 'SampHt' + else: + raise ValueError(f"Unknown instrument correction parameter key '{param_key}'.") + + value, refine_flag, min_val, max_val = correction_param + correction_changed = False + + # Set correction parameter value + if value is not None: + proj.data[hist.name]['Instrument Parameters'][0][gsas_key][1] = value + correction_changed = True + + # Set refinement flag + if refine_flag is not None: + proj.data[hist.name]['Instrument Parameters'][0][gsas_key][2] = refine_flag + + if print_info: + if correction_changed: + # Print changed parameter + print(f"Set instrument correction parameter '{param_key}' to {value} with refine_flag={refine_flag}") + + # Always print refinement flag even if value not changed + else: + print(f"No change to instrument correction parameter '{param_key}'. Refine_flag={refine_flag}") + + # Print instrument parameters set and flags if requested + if print_info: + # Check if 1st and 2nd indices are different for any instrument parameteters + # If they are all the same, then no parameters were changed + params_changed = False + for param in proj.data[hist.name]['Instrument Parameters'][0].values(): + if param[0] != param[1]: + params_changed = True + break + if not params_changed: + print("No instrument parameters were changed from their default values.") + + # print refine flags for all instrument parameters + print("Instrument parameter refinement flags:") + for param_key, param in proj.data[hist.name]['Instrument Parameters'][0].items(): + print(f" {param_key}: refine_flag={param[2]}") + + +# --- verbatim from gateways/gsasii/executors.py at 2810c82 --- + +def execute_spf_refinement( + proj: Any, + hist: Any, + recipe: RecipeModel, + verbose: bool +) -> dict: + """ + Execute single peak fitting refinement. + + This replaces the previous "peaks_only" strategy. + + Args: + proj: GSAS-II project object + hist: Histogram object + recipe: Validated recipe model + verbose: Print detailed progress + + Returns: + Result dict with success, rwp, elapsed_time + """ + controls = recipe.payload.refinement_controls + spf_mode = "useIP" if controls.single_peak_fitting_mode.use_instrument_profile else "hold" + + if verbose: + print(f"\n{'='*60}") + print(f"Executing Single Peak Fitting") + print(f" Mode: {spf_mode} ({'use instrument profile' if spf_mode == 'useIP' else 'refine peak widths'})") + print(f" Cycles: {controls.refinement_cycles}") + print(f"{'='*60}\n") + + # Execute single peak fitting (cycles already set in step 10) + peak_result = hist.refine_peaks(mode=spf_mode) + + # Extract Rwp from peak_result + rwp_final = peak_result[3].get("Rwp") if len(peak_result) > 3 else None + + if rwp_final is None: + return { + 'success': False, + 'rwp': None, + 'error': "Single peak fitting produced no Rwp — hist.refine_peaks() may have failed silently", + } + + if verbose: + print(f"Single peak fitting complete. Final Rwp: {rwp_final:.3f}%\n") + + return { + 'success': True, + 'rwp': rwp_final + } + + +# --- verbatim from gateways/gsasii/executors.py at 2810c82 --- + +SCHEMA_EXECUTORS = { + 'GSASII_Rietveld': execute_rietveld_refinement, + 'GSASII_SPF': execute_spf_refinement +} + + +# --- verbatim from gateways/gsasii/executors.py at 2810c82 --- + +def run_refinement(recipe: RecipeModel, output_dir: Path, verbose: bool = False, method: str = 'server') -> dict: + """ + Internal execution engine: run a GSAS-II refinement using already-loaded libraries. + + This function is the single execution entry point used by all paths: + - ``gsas_server.py`` calls it directly after receiving an HTTP request + (``method='server'``). + - ``GSASClient._submit_via_subprocess()`` calls it in-process as the + subprocess fallback (``method='subprocess'``). Despite the method name, + no OS subprocess is spawned; see :meth:`GSASClient._submit_via_subprocess`. + + For programmatic use, prefer ``powderline.run()`` which validates the recipe, + dispatches to the appropriate backend, and normalises all structured data + fields to pandas DataFrames. + + Args: + recipe: Validated RecipeModel instance (call ``powderline.validate()`` + first if starting from a raw dict). + output_dir: Directory where output files are written (``dummy.gpx``, + ``dummy.lst``, CSV reports, ``fit_profile.txt``). Must exist or be + creatable. + verbose: If True, print detailed progress to stdout. + method: Execution method identifier string stored in the return dict. + Callers set this to ``'server'`` or ``'subprocess'`` so downstream + consumers know how the run was executed. Default ``'server'``. + + Returns: + dict with keys: + + - ``success`` (bool) + - ``run_id`` (str) — UUID4 string + - ``rwp`` (float | None) + - ``elapsed_time`` (float) — wall-clock seconds + - ``method`` (str) — value of the ``method`` argument + - ``output_files`` (list[str]) — paths of all files written to + ``output_dir`` (informational) + - ``fit_profile`` (dict) — column-oriented dict (JSON-serializable); + normalised to ``pd.DataFrame`` by ``run()`` + - ``unit_cell_data`` (dict) — ``{phase: list-of-records}`` + - ``peak_list_data`` (dict) — ``{phase: list-of-records}`` + - ``refined_parameters`` (list) — list-of-records (9-column schema) + - ``spf_peaks`` (dict) — column-oriented dict; populated for + ``GSASII_SPF`` runs, empty dict otherwise + - ``spf_convergence_diagnostics`` (dict) — column-oriented dict; + populated when SPF peaks have convergence issues, empty dict otherwise + - ``error`` (str) — present when ``success=False`` + - ``traceback`` (str) — present when an unhandled exception occurred + + Example:: + + from powderline.kicker import run_refinement + from powderline.schema import RecipeModel + recipe = RecipeModel.model_validate(recipe_dict) + result = run_refinement(recipe, Path('output'), verbose=True) + print(f"Rwp: {result['rwp']:.3f}% [run_id={result['run_id']}]") + """ + # Assign a unique ID for this refinement run immediately, before any work begins, + # so every exit path (success, GSAS-II failure, unexpected exception) can include + # it in the returned dict and server logs. + run_id = str(uuid.uuid4()) + import time + start_time = time.time() + + try: + # 1. Get xrd_data from validated model (file-less: no sample_name) + xrd_data = recipe.payload.xrd_data + + # 2. Construct name and initialize project using standard naming + gpx_path = output_dir / OUTPUT_NAMING.gpx_filename + proj = G2.G2Project(newgpx=str(gpx_path)) + + # 3. Instrument initialization - now a list of dicts [Iparm1, Iparm2] + instrument_init = recipe.payload.instrument.initialization + + # 4. Add histogram - phase1B update: xrd_data is now a dict with arrays, instrument_init is list of dicts + try: + hist = add_powder_histogram_from_arrays( + proj=proj, + tth_array=xrd_data.tth, + intensity_array=xrd_data.Itth, + intensity_weights_array=xrd_data.Itth_weights, + histogram_name=OUTPUT_NAMING.histogram_name, + instrument_prm_dict=instrument_init[0], # Use Iparm1 (first dict) + comments=None, + phases=None, # Optionally could link to phases if order of operations + ) + except Exception as e: + import traceback as _tb + tb = _tb.format_exc() + elapsed_time = time.time() - start_time + return { + 'success': False, + 'rwp': None, + 'run_id': run_id, + 'error': f"Failed to load XRD data or instrument parameters: {str(e)}", + 'traceback': tb, + 'elapsed_time': elapsed_time, + 'method': method + } + + # After adding histogram, set hist scale and refine flag to defaults + set_hist_scale(proj, hist, hist_scale_val=DEFAULT_HIST_SCALE_VAL, hist_scale_refine_flag=DEFAULT_HIST_SCALE_REFINE_FLAG, print_info=verbose) + + # 5. Set fit range + fit_range = recipe.payload.fit_range if recipe.payload.fit_range else (None, None) + if fit_range != (None, None): + set_fit_range_hist(hist, fit_range, print_info=verbose) + + # 6. Set background - Chebyshev and single peaks + if recipe.payload.background is not None: + background_dict = recipe.payload.background.model_dump(mode='json') + if recipe.payload.background.chebyshev is not None: + chebyshev_dict = background_dict.get('chebyshev') + set_chebyshev_background(proj, hist, chebyshev_dict, print_info=verbose) + if recipe.payload.background.single_peaks is not None: + bkg_single_peaks_dict = background_dict.get('single_peaks') + set_single_peak_background(proj, hist, bkg_single_peaks_dict, print_info=verbose) + + # 6b. Set single peaks (Peak List mode for non-background peaks) + if recipe.payload.single_peaks is not None: + single_peaks_dict = recipe.payload.single_peaks.model_dump(mode='json') + # Check if dict has any actual peak data before calling setter + has_peak_data = False + for key in ['positions', 'intensities', 'pv_gaussian_sigma_sq', 'pv_lorentzian_gamma']: + if key in single_peaks_dict and single_peaks_dict[key] and len(single_peaks_dict[key]) > 0: + has_peak_data = True + break + if has_peak_data: + set_single_peaks(proj, hist, single_peaks_dict, print_info=verbose) + + # 7. Add phases (if present) + if recipe.payload.phases is not None and len(recipe.payload.phases) > 0: + phases_dict = recipe.payload.model_dump(mode='json')['phases'] + add_phases_from_dict(proj, hist, phases_dict, print_info=verbose) + + # 8. Parameterize phases + holds = set_phase_parameterization(proj, hist, phases_dict, print_info=verbose) + + # 8b. Apply Hold constraints for per-parameter refine flags + # (schema 0.26: fixed members of partially-refined DOF groups). + # The isinstance guard keeps mocked-setter unit tests inert. + if isinstance(holds, list) and holds: + proj.add_HoldConstr(holds) + if verbose: + print(f"Applied {len(holds)} Hold constraint(s): {holds}") + + # 9. Set instrument parameterization + instrument_param_dict = recipe.payload.instrument.parameterization.model_dump(mode='json') if recipe.payload.instrument.parameterization else None + if instrument_param_dict is not None: + instrument_param_changes = False + for key, value in instrument_param_dict.items(): + if isinstance(value, dict): + if any(v is not None for v in value.values()): + instrument_param_changes = True + break + elif value is not None: + instrument_param_changes = True + break + if instrument_param_changes: + set_instrument_parameterization(proj, hist, instrument_param_dict, print_info=verbose) + + # 10. Get refinement controls and set refinement cycles + controls = recipe.payload.refinement_controls + + # 10a. Validate single_peak_fitting_mode requirements (only for GSASII_SPF) + if recipe.schema_name == 'GSASII_SPF': + if controls.single_peak_fitting_mode is None: + raise ValueError( + f"Schema '{recipe.schema_name}' requires " + f"'single_peak_fitting_mode' to be configured in refinement_controls." + ) + + num_cycles = controls.refinement_cycles + if num_cycles != proj.get_Controls('cycles'): + set_refinement_cycles(proj, num_cycles, print_info=verbose) + + # 11. Execute refinement using schema-based dispatch + executor = SCHEMA_EXECUTORS.get(recipe.schema_name) + if executor is None: + raise ValueError(f"Unknown schema_name: {recipe.schema_name}. Valid options: {list(SCHEMA_EXECUTORS.keys())}") + + try: + result = executor(proj, hist, recipe, verbose) + + # Check if executor succeeded + if not result.get('success', False): + elapsed_time = time.time() - start_time + return { + 'success': False, + 'rwp': None, + 'run_id': run_id, + 'error': result.get( + 'error', + f"Refinement executor '{recipe.schema_name}' returned success=False without an error message", + ), + 'traceback': result.get('traceback'), + 'elapsed_time': elapsed_time, + 'method': method + } + + except Exception as e: + import traceback as _tb + tb = _tb.format_exc() + elapsed_time = time.time() - start_time + return { + 'success': False, + 'rwp': None, + 'run_id': run_id, + 'error': f"Refinement failed: {str(e)}", + 'traceback': tb, + 'elapsed_time': elapsed_time, + 'method': method + } + + # Save GPX file after successful refinement + proj.save(str(gpx_path)) + + # 12. Extract final Rwp from result + rwp = result.get('rwp') + + # --- Post-refinement data extraction --- + # Each helper writes its own output files and returns JSON-serializable data. + # run() normalises all dicts/lists to DataFrames for programmatic callers. + + fit_profile_data = _extract_fit_profile(hist, output_dir) + + spf_peaks_data, spf_diagnostics_data = _extract_spf_peak_report( + proj, hist, recipe, output_dir, verbose + ) + + # Extract refined parameters (needed for unit cell ESDs in _extract_phase_reports) + if verbose: + print("\nExtracting refined parameters with ESDs...") + + # NOTE: .lst fallback (extract_refined_params_from_lst) is preserved for Phase 2 + # reference but its call site is intentionally disabled: .lst parsing is incomplete + # (no atom/HAP params; different naming scheme) and obscured covariance failures. + param_dict = extract_refined_params_from_project(proj, verbose=verbose) + + # Only raise if GSAS-II populated a varyList (i.e. parameters were actually varied) + # but we still got nothing back. An empty param_dict with an empty varyList is + # expected and correct for simulation mode (refinement_cycles=1, all params locked) + # and for GSASII_SPF (no standard covariance matrix produced). + cov_vary_list = ( + proj.data.get('Covariance', {}).get('data', {}).get('varyList', []) + ) + if len(param_dict) == 0 and len(cov_vary_list) > 0: + raise RuntimeError( + "GSAS-II populated a varyList but refined parameter extraction returned " + "empty results; covariance data may be corrupt. " + "Ensure refinement completed successfully. " + "Check the .lst file for GSAS-II error details." + ) + + # --- .lst fallback (disabled) --- + # NOTE: preserved for Phase 2 reference — .lst parsing is incomplete + # (no atom/HAP params; different naming scheme); covariance extraction is the + # only supported path. + # if len(param_dict) == 0: + # lst_file = output_dir / OUTPUT_NAMING.lst_filename + # if lst_file.exists(): + # if verbose: + # print(" Covariance data not available, parsing .lst file...") + # param_dict = extract_refined_params_from_lst(lst_file) + + unit_cell_data, peak_list_data = _extract_phase_reports( + proj, hist, recipe, param_dict, output_dir + ) + + refined_parameters_data = _extract_refined_parameters( + param_dict, output_dir, proj, verbose + ) + + elapsed_time = time.time() - start_time + return { + 'success': True, + 'run_id': run_id, + 'rwp': rwp, + 'elapsed_time': elapsed_time, + 'method': method, + 'output_files': [str(f) for f in output_dir.glob('*')], + 'fit_profile': fit_profile_data, + 'unit_cell_data': unit_cell_data, + 'peak_list_data': peak_list_data, + 'refined_parameters': refined_parameters_data, + 'spf_peaks': spf_peaks_data, + 'spf_convergence_diagnostics': spf_diagnostics_data, + } + + except Exception as e: + import traceback + tb = traceback.format_exc() + import logging + logging.getLogger(__name__).error("run_refinement failed:\n%s", tb) + elapsed_time = time.time() - start_time + return { + 'success': False, + 'rwp': None, + 'run_id': run_id, + 'error': str(e), + 'traceback': tb, + 'elapsed_time': elapsed_time, + 'method': method + } + + +# --- verbatim from gateways/gsasii/extractors.py at 2810c82 --- + +def export_refined_parameters_csv( + param_dict: Dict[str, Dict[str, Any]], + output_file: Path, + proj: Any = None, + include_category: bool = True +) -> Optional[pd.DataFrame]: + """ + Export refined parameters to CSV file and return the DataFrame. + + Args: + param_dict: Dictionary from extract_refined_params_* functions + output_file: Path to output CSV file + proj: GSAS-II project object (optional, needed for phase/atom names) + include_category: If True, add 'category' and 'descriptive_name' columns (default: True) + + Returns: + DataFrame with the exported parameters, or None if param_dict is empty. + + Output CSV / DataFrame columns: + - parameter_name: GSAS-II internal parameter name + - descriptive_name: Human-readable parameter description (if include_category=True) + - phase_name: Name of associated phase (None if not phase-specific) + - phase_idx: Index of associated phase (None if not phase-specific) + - atom_name: Label of associated atom (None if not atom-specific) + - atom_idx: Index of associated atom (None if not atom-specific) + - value: Refined value + - esd: Estimated standard deviation (None for fixed parameters) + - category: Parameter category (instrument, background, cell, etc.) if include_category=True + """ + if len(param_dict) == 0: + return None + + # Build phase and atom mappings if proj is available + phase_mapping = {} + atom_mapping = {} + if proj is not None: + phase_mapping = build_phase_name_mapping(proj) + atom_mapping = build_atom_name_mapping(proj) + + rows = [] + for param_name, param_data in param_dict.items(): + # Parse parameter associations + associations = parse_parameter_associations(param_name) + phase_idx = associations["phase_idx"] + atom_idx = associations["atom_idx"] + + # Lookup names from mappings + phase_name = phase_mapping.get(phase_idx) if phase_idx is not None else None + atom_name = atom_mapping.get(phase_idx, {}).get(atom_idx) if phase_idx is not None and atom_idx is not None else None + + row = { + "parameter_name": param_name, + "value": param_data["value"], + "esd": param_data.get("esd"), + "phase_name": phase_name, + "phase_idx": phase_idx, + "atom_name": atom_name, + "atom_idx": atom_idx + } + + if include_category: + # Add descriptive name + row["descriptive_name"] = get_descriptive_param_name(param_name) + + # Determine category from parameter name + if '::A' in param_name and len(param_name.split('::')[1]) <= 2: + # Reciprocal metric tensor (A0-A5) + category = 'reciprocal_metric_tensor' + elif '::' in param_name: + # Other phase parameters (atoms, etc.) + param_type = param_name.split('::')[1] + if 'frac' in param_type: + category = 'atom_occupancy' + elif 'Uiso' in param_type: + category = 'atom_displacement_isotropic' + elif param_type.startswith('AU') and len(param_type) > 2: + category = 'atom_displacement_anisotropic' + elif param_type.startswith('A') and param_type[1] in ['x', 'y', 'z']: + category = 'atom_position' + else: + category = 'phase_other' + elif param_name.startswith('cell:'): + category = 'unit_cell' + elif ':' in param_name: + parts = param_name.split(':') + if len(parts) >= 3: + param_part = parts[2] + if param_part.startswith('Back;'): + category = 'background' + elif param_part.startswith('BkPk'): + category = 'background_peak' + elif param_part == 'Scale': + category = 'scale' + elif param_part in ['U', 'V', 'W', 'X', 'Y', 'Z']: + category = 'instrument_broadening' + elif param_part in ['Lam', 'Zero', 'SH/L', 'Polariz']: + category = 'instrument' + else: + category = 'other' + else: + category = 'other' + else: + category = 'other' + + row["category"] = category + + rows.append(row) + + # Sort by category, then parameter name + if include_category: + rows.sort(key=lambda x: (x["category"], x["parameter_name"])) + else: + rows.sort(key=lambda x: x["parameter_name"]) + + # Reorder columns: parameter_name, descriptive_name, phase_name, phase_idx, atom_name, atom_idx, value, esd, category + if include_category: + df = pd.DataFrame(rows, columns=["parameter_name", "descriptive_name", "phase_name", "phase_idx", "atom_name", "atom_idx", "value", "esd", "category"]) + else: + df = pd.DataFrame(rows, columns=["parameter_name", "phase_name", "phase_idx", "atom_name", "atom_idx", "value", "esd"]) + + # Convert index columns to nullable integer type (Int64) to avoid scientific notation formatting + df['phase_idx'] = df['phase_idx'].astype('Int64') + df['atom_idx'] = df['atom_idx'].astype('Int64') + + # Format float columns independently: values need 6 decimal places, ESDs need 8 + df['value'] = df['value'].apply(lambda x: f"{x:.6e}" if pd.notna(x) else "") + df['esd'] = df['esd'].apply(lambda x: f"{x:.8e}" if pd.notna(x) else "") + + df.to_csv(output_file, index=False, lineterminator="\n") + return df + + +# --- verbatim from gateways/gsasii/extractors.py at 2810c82 --- + +def _extract_fit_profile(hist: Any, output_dir: Path) -> dict: + """Extract fit profile arrays from histogram and save fit_profile.txt. + + Args: + hist: GSAS-II histogram object after refinement. + output_dir: Directory to write ``fit_profile.txt``. + + Returns: + dict: Column-oriented data (JSON-serializable) with keys ``two_theta``, + ``y_obs``, ``y_weights``, ``y_calc``, ``y_diff``, ``y_bkg``, ``q_values``, + ``d_spacings``. + """ + two_theta = hist.getdata(datatype="X") + q_values = hist.getdata(datatype="Q") + d_spacings = hist.getdata(datatype="d") + y_obs = hist.getdata(datatype="Yobs") + y_weights = hist.getdata(datatype="Yweight") + y_calc = hist.getdata(datatype="Ycalc") + y_bkg = hist.getdata(datatype="Background") + y_diff = hist.getdata(datatype="Residual") + + fit_profile_df = pd.DataFrame({ + "two_theta": two_theta, + "y_obs": y_obs, + "y_weights": y_weights, + "y_calc": y_calc, + "y_diff": y_diff, + "y_bkg": y_bkg, + "q_values": q_values, + "d_spacings": d_spacings, + }) + fit_profile_df.to_csv( + output_dir / "fit_profile.txt", sep="\t", float_format="%.8f", + header=True, index=False, lineterminator="\n" + ) + return {col: fit_profile_df[col].tolist() for col in fit_profile_df.columns} + + +# --- verbatim from gateways/gsasii/extractors.py at 2810c82 --- + +def _extract_spf_peak_report( + proj: Any, hist: Any, recipe: RecipeModel, output_dir: Path, verbose: bool +) -> tuple[dict, dict]: + """Extract single peak fitting results and save report files. + + Called only for ``GSASII_SPF`` runs where ``recipe.payload.single_peaks`` + is set. Returns two column-oriented dicts (JSON-serializable) that are + normalised to DataFrames by ``run()``. + + Also writes: + - ``single_peaks_report.txt`` — per-peak widths and convergence status + - ``peak_convergence_diagnostics.txt`` — only when peaks have issues + + Args: + proj: GSAS-II project object after refinement. + hist: GSAS-II histogram object. + recipe: Validated RecipeModel. + output_dir: Directory to write report files. + verbose: If True, print convergence warnings to stdout. + + Returns: + ``(spf_peaks_data, spf_diagnostics_data)`` where each is a + column-oriented dict, or ``({}, {})`` when single peaks are not used. + """ + if recipe.payload.single_peaks is None: + return {}, {} + + peak_list = proj.data.get(hist.name, {}).get('Peak List', {}).get('peaks', []) + if not peak_list: + return {}, {} + + peak_data = [] + convergence_diagnostics = [] + converged_count = 0 + aphysical_count = 0 + + for i, peak in enumerate(peak_list): + pos, pos_flag, intensity, int_flag, sig_sq, sig_sq_flag, gamma, gam_flag = peak + + # Determine convergence status + status = "converged" + if np.isnan(sig_sq) or np.isnan(gamma): + status = "NaN_failed" + elif sig_sq <= 0 and gamma < 0: + status = "negative_sigma_sq_and_gamma_warning" + aphysical_count += 1 + elif sig_sq <= 0: + status = "zero_or_negative_sigma_sq_warning" + aphysical_count += 1 + elif gamma < 0: + status = "negative_gamma_warning" + aphysical_count += 1 + else: + converged_count += 1 + + # Copy GSAS-II behavior for PV calc — use a small positive default for + # aphysical (zero/negative) sigma_sq or gamma values + sigma = np.sqrt(sig_sq) if sig_sq > 0 else DEFAULT_SPF_SIGMA_MIN + # NOTE: GSAS-II gamma = FWHM_L (not HWHM), so we use gamma directly + gamma_calc = gamma if gamma > 0 else DEFAULT_SPF_GAMMA_MIN + + # Calculate widths in degrees (GSAS-II stores in centidegrees) + fwhm_g, fwhm_l, fwhm_pv, ib_g, ib_l, ib_pv, valid, warning_msg = calculate_peak_widths( + sigma / 100, gamma_calc * 0.5 / 100 + ) + + # Get GSAS-II verification (may fail for aphysical values) + try: + # GSAS-II uses gamma directly as FWHM_L (differs from scipy, TOPAS, literature) + fwhm_gsas = G2pwd.getgamFW(gamma_calc / 100, sigma / 100) + except Exception: + fwhm_gsas = np.nan + + converged_bool = status == "converged" + peak_data.append([ + pos, intensity, sigma, sig_sq, gamma_calc, + fwhm_g, fwhm_l, fwhm_pv, ib_g, ib_l, ib_pv, + fwhm_gsas, converged_bool, status + ]) + + if status != "converged": + convergence_diagnostics.append({ + 'peak_index': i, + 'position_2theta': pos, + 'final_sigma_sq': sig_sq, + 'final_gamma': gamma, + 'status': status, + 'notes': warning_msg if warning_msg else "" + }) + + # Build DataFrames + spf_cols = [ + "position_2theta", "intensity", "sigma", "sigma_squared", "gamma", + "fwhm_gaussian", "fwhm_lorentzian", "fwhm_pseudovoigt", + "integral_breadth_gaussian", "integral_breadth_lorentzian", + "integral_breadth_pseudovoigt", "fwhm_gsas_verification", + "converged", "convergence_detail" + ] + spf_peaks_df = pd.DataFrame(peak_data, columns=spf_cols) + diag_df = pd.DataFrame(convergence_diagnostics) if convergence_diagnostics else pd.DataFrame( + columns=['peak_index', 'position_2theta', 'final_sigma_sq', 'final_gamma', 'status', 'notes'] + ) + + # Write single_peaks_report.txt (custom mixed-type formatting) + header_comment = ( + f"{converged_count} of {len(peak_list)} peaks converged; " + f"{aphysical_count} peaks have aphysical values (warnings)" + ) + header_cols = "\t".join(spf_cols) + peak_array = np.array(peak_data, dtype=object) + with open(output_dir / "single_peaks_report.txt", 'w', newline="\n") as f: + f.write(f"# {header_comment}\n") + f.write(f"{header_cols}\n") + for row in peak_array: + formatted_row = [] + for j, val in enumerate(row): + col_name = spf_cols[j] + if col_name == "converged": + formatted_row.append(str(val).lower()) + elif col_name == "convergence_detail": + formatted_row.append(str(val)) + else: + formatted_row.append("nan" if np.isnan(val) else f"{val:.8f}") + f.write("\t".join(formatted_row) + "\n") + + # Write diagnostics file if there are issues + if convergence_diagnostics: + diag_df.to_csv( + output_dir / "peak_convergence_diagnostics.txt", sep="\t", + float_format="%.8f", header=True, index=False, lineterminator="\n" + ) + if verbose: + print( + f"\nWarning: {len(convergence_diagnostics)} of {len(peak_list)} " + f"peaks have convergence issues" + ) + print(f" See {output_dir / 'peak_convergence_diagnostics.txt'} for details\n") + + return ( + {col: spf_peaks_df[col].tolist() for col in spf_peaks_df.columns}, + {col: diag_df[col].tolist() for col in diag_df.columns}, + ) + + +# --- verbatim from gateways/gsasii/extractors.py at 2810c82 --- + +def _extract_phase_reports( + proj: Any, hist: Any, recipe: RecipeModel, param_dict: dict, output_dir: Path +) -> tuple[dict, dict]: + """Extract unit cell and peak list reports for all phases. + + Writes ``{phase}_unit_cell_report.csv`` and ``{phase}_peak_list_report.csv`` + for each phase. + + Args: + proj: GSAS-II project object after refinement. + hist: GSAS-II histogram object. + recipe: Validated RecipeModel (used to check if phases exist). + param_dict: Refined parameter dict from + ``extract_refined_params_from_project()`` (needed for cell ESDs). + output_dir: Directory to write CSV files. + + Returns: + ``(unit_cell_data, peak_list_data)`` — each is a + ``{phase_name: list-of-records}`` dict (JSON-serializable). + """ + unit_cell_data: dict = {} + peak_list_data: dict = {} + + if recipe.payload.phases is None or len(recipe.payload.phases) == 0: + return unit_cell_data, peak_list_data + + phase_names = [p.name for p in proj.phases() if p is not None] + + for phase_idx, phase_name in enumerate(phase_names): + try: + unit_cell = proj.data['Phases'][phase_name]['General']['Cell'][1:8] + except (KeyError, IndexError): + continue + + cell_params = ["cell_a", "cell_b", "cell_c", "cell_alpha", "cell_beta", "cell_gamma", "cell_volume"] + esds = calculate_cell_esds_from_A_matrix(phase_idx, proj, phase_name) + + # NOTE: .lst fallback for cell ESDs is disabled — preserved for Phase 2 reference. + # calculate_cell_esds_from_A_matrix raises RuntimeError on failure rather than + # returning [None]*7, so the fallback block below is unreachable. + # if all(esd is None for esd in esds): + # cell_param_keys = ['a', 'b', 'c', 'alpha', 'beta', 'gamma', 'volume'] + # esds = [param_dict.get(f"cell:{k}", {}).get("esd") for k in cell_param_keys] + + unit_cell_df = pd.DataFrame({ + "parameter": cell_params, + "value": unit_cell, + "esd": esds, + }) + unit_cell_df.to_csv( + output_dir / f"{phase_name}_unit_cell_report.csv", + float_format="%.8f", index=False, lineterminator="\n" + ) + unit_cell_data[phase_name] = json.loads(unit_cell_df.to_json(orient='records')) + + # Peak list for each phase + reflection_lists = proj.data.get(hist.name, {}).get('Reflection Lists', {}) + for phase_name in phase_names: + if phase_name not in reflection_lists: + continue + phase_ref = reflection_lists[phase_name] + if 'RefList' not in phase_ref: + continue + reflection_list = phase_ref['RefList'] + headers = [ + "h", "k", "l", "multiplicity", "d_spacing", "2theta", + "sigma_squared", "gamma", "F_obs_squared", "F_calc_squared", + "phase", "I_corr", "Prfo", "Trans", "ExtP" + ] + peak_list_df = pd.DataFrame(reflection_list, columns=headers) + peak_list_df.to_csv( + output_dir / f"{phase_name}_peak_list_report.csv", + float_format="%.8f", index=False, lineterminator="\n" + ) + peak_list_data[phase_name] = json.loads(peak_list_df.to_json(orient='records')) + + return unit_cell_data, peak_list_data + + +# --- verbatim from gateways/gsasii/extractors.py at 2810c82 --- + +def _extract_refined_parameters( + param_dict: dict, output_dir: Path, proj: Any, verbose: bool +) -> list: + """Export refined parameters to CSV and return as list-of-records. + + Args: + param_dict: From ``extract_refined_params_from_project()``. + output_dir: Directory to write ``refined_parameters.csv``. + proj: GSAS-II project object (for phase/atom name mappings). + verbose: If True, print export status to stdout. + + Returns: + List of dicts (records) with the 9-column schema. Returns ``[]`` + when ``param_dict`` is empty (simulation mode, SPF, etc.). + """ + if not param_dict: + if verbose: + print(" No refined parameters found to export") + return [] + + try: + refined_params_csv = output_dir / "refined_parameters.csv" + refined_params_df = export_refined_parameters_csv( + param_dict, refined_params_csv, proj=proj + ) + if refined_params_df is not None: + if verbose: + print(f" Exported {len(param_dict)} refined parameters to {refined_params_csv.name}") + return refined_params_df.to_dict('records') + except Exception as e: + if verbose: + print(f" Warning: Failed to export refined parameters CSV: {e}") + return [] diff --git a/src/powderline/gateways/gsasii/client.py b/src/powderline/gateways/gsasii/client.py index 5b4bd86..0cefaf4 100644 --- a/src/powderline/gateways/gsasii/client.py +++ b/src/powderline/gateways/gsasii/client.py @@ -154,6 +154,12 @@ def submit_simulation(self, recipe, output_dir: Path, False, returns a structured error explaining the filesystem-view mismatch. """ + result = self._dispatch(recipe, output_dir, verbose, auto_start_server, use_server) + return _with_validation_warnings(recipe, result) + + def _dispatch(self, recipe, output_dir: Path, verbose: bool, auto_start_server: bool, + use_server: bool) -> Dict[str, Any]: + """Run the recipe on the server or in this process (see :meth:`submit_simulation`).""" recipe_dict = self._normalize_recipe(recipe) if not use_server: @@ -235,7 +241,8 @@ def submit_simulation(self, recipe, output_dir: Path, 'method': 'none' } - def _normalize_recipe(self, recipe) -> dict: + @staticmethod + def _normalize_recipe(recipe) -> dict: """Normalize recipe input to a plain JSON-serializable dict. Accepts Path (loaded from disk), dict (used as-is), or RecipeModel @@ -385,13 +392,13 @@ def _submit_via_subprocess(self, recipe_dict: dict, output_dir: Path, """ import time from powderline.kicker import run_refinement - from powderline.schema import RecipeModel + from powderline.gateways.gsasii.routing import model_from_request from pydantic import ValidationError output_dir.mkdir(parents=True, exist_ok=True) try: - recipe = RecipeModel.model_validate(recipe_dict) + recipe = model_from_request(recipe_dict) # native gsasii.* or 0.26.0 (re/04) except ValidationError as e: error_lines = [ f"{' -> '.join(str(loc) for loc in err['loc'])}: {err['msg']}" @@ -452,3 +459,31 @@ def run_simulation(recipe, output_dir: Path, """ client = GSASClient() return client.submit_simulation(recipe, output_dir, verbose, auto_start_server) + + +def _with_validation_warnings(recipe, result: Dict[str, Any]) -> Dict[str, Any]: + """Prepend a native ``gsasii.*`` recipe's validation warnings to the run's ``warnings`` (A115). + + They are collected from the recipe as given (the model the caller validated, + or the dict / file validated here), so a default the recipe left out is + reported even though the server receives the dumped model, which states it. + A 0.26.0 recipe, or one that does not validate (the run already says why), + gets none. + """ + from pydantic import ValidationError + from powderline.gateways.gsasii.schema import is_native_recipe, validate_recipe + from powderline.schema_core import collect_warnings + + if isinstance(recipe, Path): + try: + recipe = GSASClient._normalize_recipe(recipe) + except (OSError, ValueError): + return result + if not is_native_recipe(recipe): + return result + try: + found = collect_warnings(validate_recipe(recipe)) + except ValidationError: + return result + result['warnings'] = found + list(result.get('warnings') or []) + return result diff --git a/src/powderline/gateways/gsasii/executors.py b/src/powderline/gateways/gsasii/executors.py index 60acfe2..03fcd75 100644 --- a/src/powderline/gateways/gsasii/executors.py +++ b/src/powderline/gateways/gsasii/executors.py @@ -25,6 +25,10 @@ from powderline.gateways.gsasii.constraints import atom_refinement_plan, cell_refinement_plan from powderline._status import CHECK, CROSS, INFO, WARN, emoji from dataclasses import dataclass +from powderline.gateways.gsasii.schema import GsasiiRietveldRecipe, GsasiiSpfRecipe, gsasii_iparm1 +from powderline.gateways.gsasii.fit_report import fit_report +from powderline.gateways.gsasii.sites import check_sites +from powderline.exceptions import EngineExecutionError from powderline.gateways.gsasii.helpers import ( DEFAULT_HIST_SCALE_REFINE_FLAG, @@ -177,7 +181,15 @@ def execute_spf_refinement( print(f"{'='*60}\n") # Execute single peak fitting (cycles already set in step 10) - peak_result = hist.refine_peaks(mode=spf_mode) + try: + peak_result = hist.refine_peaks(mode=spf_mode) + except TypeError as exc: # DoPeakFit failed and returned None; refine_peaks then indexes it (EB-38) + if "'NoneType' object is not subscriptable" not in str(exc): # any other TypeError keeps its own traceback + raise + raise EngineExecutionError( + "GSAS-II's peak fit failed: DoPeakFit returned no result (GSAS-II prints the reason only in its " + "debug mode; GSAS-II quirk EB-38). Check the peak positions, widths and fit range." + ) from exc # Extract Rwp from peak_result rwp_final = peak_result[3].get("Rwp") if len(peak_result) > 3 else None @@ -194,14 +206,17 @@ def execute_spf_refinement( return { 'success': True, - 'rwp': rwp_final + 'rwp': rwp_final, + 'engine_rvals': peak_result[3], # DoPeakFit's own residuals, for engine_details (A23) } # Schema executor registry SCHEMA_EXECUTORS = { 'GSASII_Rietveld': execute_rietveld_refinement, - 'GSASII_SPF': execute_spf_refinement + 'GSASII_SPF': execute_spf_refinement, + 'gsasii.rietveld': execute_rietveld_refinement, + 'gsasii.spf': execute_spf_refinement, } @@ -209,7 +224,7 @@ def execute_spf_refinement( # Kicker script for PowderLine ######################################## -def run_refinement(recipe: RecipeModel, output_dir: Path, verbose: bool = False, method: str = 'server') -> dict: +def run_refinement(recipe: GsasiiRietveldRecipe | GsasiiSpfRecipe, output_dir: Path, verbose: bool = False, method: str = 'server') -> dict: """ Internal execution engine: run a GSAS-II refinement using already-loaded libraries. @@ -225,8 +240,9 @@ def run_refinement(recipe: RecipeModel, output_dir: Path, verbose: bool = False, fields to pandas DataFrames. Args: - recipe: Validated RecipeModel instance (call ``powderline.validate()`` - first if starting from a raw dict). + recipe: Validated ``gsasii.rietveld`` / ``gsasii.spf`` recipe model (call + ``powderline.validate()`` first if starting from a raw dict). 0.26.0 + ``RecipeModel`` recipes run through ``_legacy_executor.run_refinement``. output_dir: Directory where output files are written (``dummy.gpx``, ``dummy.lst``, CSV reports, ``fit_profile.txt``). Must exist or be creatable. @@ -259,9 +275,8 @@ def run_refinement(recipe: RecipeModel, output_dir: Path, verbose: bool = False, Example:: - from powderline.kicker import run_refinement - from powderline.schema import RecipeModel - recipe = RecipeModel.model_validate(recipe_dict) + from powderline.kicker import run_refinement, validate + recipe = validate(recipe_dict) # a gsasii.rietveld recipe dict result = run_refinement(recipe, Path('output'), verbose=True) print(f"Rwp: {result['rwp']:.3f}% [run_id={result['run_id']}]") """ @@ -280,8 +295,8 @@ def run_refinement(recipe: RecipeModel, output_dir: Path, verbose: bool = False, gpx_path = output_dir / OUTPUT_NAMING.gpx_filename proj = G2.G2Project(newgpx=str(gpx_path)) - # 3. Instrument initialization - now a list of dicts [Iparm1, Iparm2] - instrument_init = recipe.payload.instrument.initialization + # 3. Instrument parameters: GSAS-II's Iparm1 built from the instrument block (A99, A103) + iparm1 = gsasii_iparm1(recipe.payload.instrument) # 4. Add histogram - phase1B update: xrd_data is now a dict with arrays, instrument_init is list of dicts try: @@ -291,7 +306,7 @@ def run_refinement(recipe: RecipeModel, output_dir: Path, verbose: bool = False, intensity_array=xrd_data.Itth, intensity_weights_array=xrd_data.Itth_weights, histogram_name=OUTPUT_NAMING.histogram_name, - instrument_prm_dict=instrument_init[0], # Use Iparm1 (first dict) + instrument_prm_dict=iparm1, comments=None, phases=None, # Optionally could link to phases if order of operations ) @@ -313,7 +328,7 @@ def run_refinement(recipe: RecipeModel, output_dir: Path, verbose: bool = False, set_hist_scale(proj, hist, hist_scale_val=DEFAULT_HIST_SCALE_VAL, hist_scale_refine_flag=DEFAULT_HIST_SCALE_REFINE_FLAG, print_info=verbose) # 5. Set fit range - fit_range = recipe.payload.fit_range if recipe.payload.fit_range else (None, None) + fit_range = (recipe.payload.fit_range.min, recipe.payload.fit_range.max) if recipe.payload.fit_range else (None, None) if fit_range != (None, None): set_fit_range_hist(hist, fit_range, print_info=verbose) @@ -328,7 +343,7 @@ def run_refinement(recipe: RecipeModel, output_dir: Path, verbose: bool = False, set_single_peak_background(proj, hist, bkg_single_peaks_dict, print_info=verbose) # 6b. Set single peaks (Peak List mode for non-background peaks) - if recipe.payload.single_peaks is not None: + if getattr(recipe.payload, 'single_peaks', None) is not None: # gsasii.spf only single_peaks_dict = recipe.payload.single_peaks.model_dump(mode='json') # Check if dict has any actual peak data before calling setter has_peak_data = False @@ -340,9 +355,10 @@ def run_refinement(recipe: RecipeModel, output_dir: Path, verbose: bool = False, set_single_peaks(proj, hist, single_peaks_dict, print_info=verbose) # 7. Add phases (if present) - if recipe.payload.phases is not None and len(recipe.payload.phases) > 0: + if getattr(recipe.payload, 'phases', None) is not None and len(recipe.payload.phases) > 0: # gsasii.rietveld only phases_dict = recipe.payload.model_dump(mode='json')['phases'] add_phases_from_dict(proj, hist, phases_dict, print_info=verbose) + check_sites(proj, recipe) # GSAS-II's reading of each atom's site (A117; EB-40, EB-05) # 8. Parameterize phases holds = set_phase_parameterization(proj, hist, phases_dict, print_info=verbose) @@ -356,7 +372,7 @@ def run_refinement(recipe: RecipeModel, output_dir: Path, verbose: bool = False, print(f"Applied {len(holds)} Hold constraint(s): {holds}") # 9. Set instrument parameterization - instrument_param_dict = recipe.payload.instrument.parameterization.model_dump(mode='json') if recipe.payload.instrument.parameterization else None + instrument_param_dict = recipe.payload.instrument.model_dump(mode='json') # single source (A99, A103) if instrument_param_dict is not None: instrument_param_changes = False for key, value in instrument_param_dict.items(): @@ -374,7 +390,7 @@ def run_refinement(recipe: RecipeModel, output_dir: Path, verbose: bool = False, controls = recipe.payload.refinement_controls # 10a. Validate single_peak_fitting_mode requirements (only for GSASII_SPF) - if recipe.schema_name == 'GSASII_SPF': + if recipe.schema_name == 'gsasii.spf': if controls.single_peak_fitting_mode is None: raise ValueError( f"Schema '{recipe.schema_name}' requires " @@ -496,6 +512,9 @@ def run_refinement(recipe: RecipeModel, output_dir: Path, verbose: bool = False, 'refined_parameters': refined_parameters_data, 'spf_peaks': spf_peaks_data, 'spf_convergence_diagnostics': spf_diagnostics_data, + # Standard fit statistics from core fitstats; 'rwp' becomes the standard Rwp and + # GSAS-II's own values go to engine_details (A23, A45, A125) + **fit_report(proj, hist, recipe, engine_rwp=rwp, engine_rvals=result.get('engine_rvals')), } except Exception as e: diff --git a/src/powderline/gateways/gsasii/extractors.py b/src/powderline/gateways/gsasii/extractors.py index 737e27b..a3862c9 100644 --- a/src/powderline/gateways/gsasii/extractors.py +++ b/src/powderline/gateways/gsasii/extractors.py @@ -26,6 +26,13 @@ from powderline._status import CHECK, CROSS, INFO, WARN, emoji from dataclasses import dataclass +from powderline.gateways.gsasii.report_rows import pandas_csv_rows +from powderline.reports import ( + write_fit_profile, + write_peak_list_report, + write_refined_parameters, + write_unit_cell_report, +) from powderline.gateways.gsasii.helpers import ( DEFAULT_SPF_GAMMA_MIN, DEFAULT_SPF_SIGMA_MIN, @@ -706,7 +713,10 @@ def export_refined_parameters_csv( df['value'] = df['value'].apply(lambda x: f"{x:.6e}" if pd.notna(x) else "") df['esd'] = df['esd'].apply(lambda x: f"{x:.8e}" if pd.notna(x) else "") - df.to_csv(output_file, index=False, lineterminator="\n") + if include_category and Path(output_file).name == "refined_parameters.csv": # the core writer's file (A83) + write_refined_parameters(Path(output_file).parent, pandas_csv_rows(df)) + else: + df.to_csv(output_file, index=False, lineterminator="\n") return df @@ -745,10 +755,7 @@ def _extract_fit_profile(hist: Any, output_dir: Path) -> dict: "q_values": q_values, "d_spacings": d_spacings, }) - fit_profile_df.to_csv( - output_dir / "fit_profile.txt", sep="\t", float_format="%.8f", - header=True, index=False, lineterminator="\n" - ) + write_fit_profile(output_dir, pandas_csv_rows(fit_profile_df, "%.8f")) # core writer (A83) return {col: fit_profile_df[col].tolist() for col in fit_profile_df.columns} @@ -776,7 +783,7 @@ def _extract_spf_peak_report( ``(spf_peaks_data, spf_diagnostics_data)`` where each is a column-oriented dict, or ``({}, {})`` when single peaks are not used. """ - if recipe.payload.single_peaks is None: + if getattr(recipe.payload, 'single_peaks', None) is None: # gsasii.spf only return {}, {} peak_list = proj.data.get(hist.name, {}).get('Peak List', {}).get('peaks', []) @@ -919,7 +926,7 @@ def _extract_phase_reports( unit_cell_data: dict = {} peak_list_data: dict = {} - if recipe.payload.phases is None or len(recipe.payload.phases) == 0: + if getattr(recipe.payload, 'phases', None) is None or len(recipe.payload.phases) == 0: # gsasii.rietveld only return unit_cell_data, peak_list_data phase_names = [p.name for p in proj.phases() if p is not None] @@ -945,10 +952,7 @@ def _extract_phase_reports( "value": unit_cell, "esd": esds, }) - unit_cell_df.to_csv( - output_dir / f"{phase_name}_unit_cell_report.csv", - float_format="%.8f", index=False, lineterminator="\n" - ) + write_unit_cell_report(output_dir, phase_name, pandas_csv_rows(unit_cell_df, "%.8f")) # core writer (A83) unit_cell_data[phase_name] = json.loads(unit_cell_df.to_json(orient='records')) # Peak list for each phase @@ -966,10 +970,8 @@ def _extract_phase_reports( "phase", "I_corr", "Prfo", "Trans", "ExtP" ] peak_list_df = pd.DataFrame(reflection_list, columns=headers) - peak_list_df.to_csv( - output_dir / f"{phase_name}_peak_list_report.csv", - float_format="%.8f", index=False, lineterminator="\n" - ) + write_peak_list_report(output_dir, phase_name, pandas_csv_rows(peak_list_df, "%.8f"), + columns=headers) # core writer (A83) peak_list_data[phase_name] = json.loads(peak_list_df.to_json(orient='records')) return unit_cell_data, peak_list_data diff --git a/src/powderline/gateways/gsasii/fit_report.py b/src/powderline/gateways/gsasii/fit_report.py new file mode 100644 index 0000000..14ecba5 --- /dev/null +++ b/src/powderline/gateways/gsasii/fit_report.py @@ -0,0 +1,150 @@ +"""GSAS-II gateway — standard fit statistics, engine details and engine warnings after a run (A23, A41, A45, A125). + +Read after the refinement from what GSAS-II stored: + +- the profile arrays (``hist.data['data'][1]``: 2theta, yobs, weights, ycalc) + and the fit window (``Limits``; excluded regions are masked), for core + :mod:`powderline.fitstats` (standard keys ``rwp``, ``r_exp``, ``gof``, + ``chi2_red``); +- ``P``, the varied-parameter count GSAS-II used: Rietveld + ``Covariance['data']['Rvals']['Nvars']`` (it includes parameters GSAS-II then + dropped as singular, EB-33), Peak List fit ``len(Peak List['sigDict'])``; +- GSAS-II's own values (``Rwp``; ``GOF``, which for a Rietveld refinement is + sqrt(reduced chi^2) and for a peak fit the reduced chi^2 itself, EB-52; + ``chisq``; ``Nobs``) + under ``engine_details``, with ``parameters_requested`` (the recipe's + independently refined parameters) next to ``parameters_varied`` (P); +- GSAS-II's refinement message (``Rvals['msg']``: parameters dropped as + singular, SVD problems, limits), which ``G2strMain.Refine`` only prints + (EB-32), as a structured warning; and a warning when GSAS-II varied fewer + parameters than the recipe requested (A128). + +A calculated pattern that is not finite inside the fit window means the +refinement diverged; GSAS-II does not report that as a failure (its wR shows +100 %, EB-52), so it is raised here as an ``EngineExecutionError`` (A130). A +simulation (``refinement_cycles == 1``) gets the same statistics as a +refinement, and every result says which it was (``simulation_mode``, A129). +A statistic that is undefined (no degrees of freedom) is ``None``. + +Runtime layer (called with GSAS-II objects); imports no GSAS-II module itself. +""" + +from __future__ import annotations + +import re + +import numpy as np +from pydantic import BaseModel + +from powderline.exceptions import EngineExecutionError, StructuredWarning +from powderline.fitstats import compute_fit_statistics +from powderline.schema_core import ChebyshevBackground, Phase, RefinableParameter +from powderline.symmetry import cell_tie_groups, coupling_groups + + +def parameters_requested(recipe) -> int: + """Independently refined parameters the recipe asks for: one per symmetry tie group (A94), per flag otherwise.""" + count = 0 + + def walk(obj) -> None: + nonlocal count + if isinstance(obj, RefinableParameter): + count += obj.refine_flag + elif isinstance(obj, ChebyshevBackground): + count += obj.num_coefficients if obj.refine_flag else 0 + elif isinstance(obj, dict): + for v in obj.values(): + walk(v) + elif isinstance(obj, list): + for v in obj: + walk(v) + elif isinstance(obj, BaseModel): + for name in type(obj).model_fields: + walk(getattr(obj, name)) + if isinstance(obj, Phase): # a refined tie group is one engine parameter + tied = [(obj.unit_cell, g.members) for g in cell_tie_groups(obj.space_group).groups] + for atom in obj.atoms.values(): + ties = coupling_groups(obj.space_group, (atom.x.value, atom.y.value, atom.z.value)) + tied += [(atom, g.members) for g in ties.xyz.groups] + if atom.Uaniso is not None: + tied += [(atom.Uaniso, g.members) for g in ties.uij.groups] + for holder, members in tied: + refined = sum(getattr(holder, m).refine_flag for m in members) + count -= max(refined - 1, 0) + + walk(recipe.payload) + return count + + +def _dropped(msg: str) -> int | None: + match = re.search(r"(\d+)\s+Parameter\(s\) dropped", msg or "") + return int(match.group(1)) if match else None + + +def _finite_or_none(v) -> float | None: + """A plain float (the result crosses JSON in server mode), or ``None`` if undefined (NaN) or absent.""" + return None if v is None or not np.isfinite(v) else float(v) + + +def fit_report(proj, hist, recipe, engine_rwp, engine_rvals=None) -> dict: + """``rwp``, ``r_exp``, ``gof``, ``chi2_red``, ``simulation_mode``, ``engine_details`` and ``warnings``. + + ``engine_rwp`` is the Rwp GSAS-II reported for the run (the executor's + value); ``engine_rvals`` the peak fit's own residuals (``DoPeakFit``'s + ``Rvals``; single peak fitting only). Raises ``EngineExecutionError`` when + the refinement diverged (a non-finite calculated pattern in the fit window). + """ + x, yobs, weights, ycalc = (hist.data["data"][1][i] for i in range(4)) + lo, hi = hist.data["Limits"][1] + xs = np.ma.getdata(x) + mask = (xs >= lo) & (xs <= hi) & ~np.ma.getmaskarray(x) + yc = np.ma.getdata(ycalc) + diverged = int(np.count_nonzero(~np.isfinite(yc[mask]))) + if diverged: + raise EngineExecutionError( + f"GSAS-II's refinement diverged: the calculated pattern is not finite at {diverged} of " + f"{int(np.count_nonzero(mask))} points in the fit window, so no result is reported (GSAS-II itself " + "does not flag this; its wR shows 100 %, EB-52). Refine fewer or less correlated parameters, or start " + "closer to the solution; the GSAS-II files in the output directory show the diverged parameters.") + spf = recipe.schema_name == "gsasii.spf" + if spf: + rvals = {k: v for k, v in (engine_rvals or {}).items() if k in ("GOF",)} + n_params = len(hist.data.get("Peak List", {}).get("sigDict") or {}) + else: + rvals = proj.data.get("Covariance", {}).get("data", {}).get("Rvals", {}) or {} + n_params = int(rvals.get("Nvars", 0)) + stats = compute_fit_statistics(np.ma.getdata(yobs), yc, np.ma.getdata(weights), n_params, mask=mask) + message = " ".join(str(rvals.get("msg", "")).split()) + num = _finite_or_none + requested = parameters_requested(recipe) + details = { + "engine": "GSAS-II", + "rwp": num(engine_rwp), + "gof": num(rvals.get("GOF")), + "chi2": num(rvals.get("chisq")), + "n_obs": None if rvals.get("Nobs") is None else int(rvals["Nobs"]), + "n_points": stats.n_points, + "parameters_requested": requested, + "parameters_varied": n_params, + "message": message or None, + } + warnings = [] + if message: + dropped = _dropped(message) + warnings.append(StructuredWarning( + code="gsasii_refinement_message", + message=(f"GSAS-II reported: {message}" + + (f" ({dropped} parameter(s) dropped as singular are still counted in " + "parameters_varied, GSAS-II's own count, EB-33)" if dropped else "")), + field_path=None)) + if n_params < requested: + warnings.append(StructuredWarning( + code="gsasii_parameters_not_varied", + message=(f"GSAS-II varied {n_params} parameter(s), but the recipe refines {requested} (one per " + "symmetry tie group); GSAS-II left the others fixed or dropped them" + + (" as singular without reporting which (single peak fitting, EB-53)" if spf else "") + + "; compare refined_parameters with the recipe's refine flags"), + field_path=None)) + return {"rwp": num(stats.rwp), "r_exp": num(stats.rexp), "gof": num(stats.gof), "chi2_red": num(stats.chi2_red), + "simulation_mode": recipe.payload.refinement_controls.refinement_cycles == 1, + "engine_details": details, "warnings": warnings} diff --git a/src/powderline/gateways/gsasii/gateway.py b/src/powderline/gateways/gsasii/gateway.py index 5ea97a1..b4367f4 100644 --- a/src/powderline/gateways/gsasii/gateway.py +++ b/src/powderline/gateways/gsasii/gateway.py @@ -1,7 +1,9 @@ """Engine-free factory for the GSAS-II gateway (entry point ``powderline.gateways:gsasii``). - **Schema layer** — ``validate`` / ``capabilities``: engine-free; runs with - GSAS-II not installed (``gateways.gsasii.validation`` is recipe-only). + GSAS-II not installed. ``gsasii.*`` recipes validate against the native + schemas (``gateways.gsasii.schema``, A39); 0.26.0 ``GSASII_*`` recipes keep + the legacy validation (``gateways.gsasii.validation``) until re/07. - **Runtime layer** — ``run``: a thin wrapper over ``powderline.kicker.run`` (behavior unchanged), imported inside the function. On a real run it first reads the installed GSAS-II version and enforces ``ENGINE_VERSION_SPEC``. @@ -14,25 +16,31 @@ from __future__ import annotations from powderline.exceptions import EngineNotAvailableError +from powderline.gateways.gsasii import schema from powderline.gateways.gsasii.validation import validate as _validate_recipe from powderline.registry import Gateway, check_engine_version NAME = "gsasii" #: Supported GSAS-II runtime versions (A41/A44); pixi pins align in re/08. ENGINE_VERSION_SPEC = "==5.7.9" -#: Workflows served (0.26.0 schema names until re/04 introduces ``gsasii.*``). -WORKFLOWS = ("GSASII_Rietveld", "GSASII_SPF") +#: Workflows served: the native schemas, and the 0.26.0 names until re/07 removes them. +WORKFLOWS = (*schema.SCHEMA_NAMES, "GSASII_Rietveld", "GSASII_SPF") #: pixi environment providing GSAS-II + how to install it (names final in re/08). ENV = "default" INSTALL_COMMAND = "pixi install" def capabilities() -> dict: - return {"name": NAME, "workflows": WORKFLOWS, "engine_version_spec": ENGINE_VERSION_SPEC} + return {"name": NAME, "workflows": WORKFLOWS, "engine_version_spec": ENGINE_VERSION_SPEC, + **schema.DECLARATIONS} def validate(recipe, *, verbose: bool = False): - """Validate a recipe (dict or ``RecipeModel``); return the ``RecipeModel`` (A50, A59).""" + """Validate a recipe (dict or model); return the validated model (A50, A59). + + A ``gsasii.*`` recipe gives its native model (``GsasiiRietveldRecipe`` / + ``GsasiiSpfRecipe``); a 0.26.0 one the legacy ``RecipeModel``. + """ return _validate_recipe(recipe, verbose=verbose) diff --git a/src/powderline/gateways/gsasii/helpers.py b/src/powderline/gateways/gsasii/helpers.py index 5eb37f6..ab2a27d 100644 --- a/src/powderline/gateways/gsasii/helpers.py +++ b/src/powderline/gateways/gsasii/helpers.py @@ -51,14 +51,16 @@ class OutputNamingConfig: def unpack_refinement_parameter(param: list | None, param_name: str = "parameter") -> tuple[Any, Any, Any, Any]: """ - Unpack refinement parameter from [value, refine_flag, min, max] format. + Unpack a refinement parameter: [value, refine_flag, min, max] (schema 0.26.0) or + [value, refine_flag] (gsasii 1.0.0, no bounds: min and max are returned as None). This utility reduces boilerplate for the common pattern of unpacking refinement parameters throughout the codebase. Use only in straightforward cases - complex conditional logic should inline the unpacking for clarity. Args: - param: RefinementParameter list [value, refine_flag, min, max] or None + param: RefinementParameter list [value, refine_flag, min, max] or + [value, refine_flag], or None param_name: Name of parameter for error messages (optional) Returns: @@ -75,10 +77,13 @@ def unpack_refinement_parameter(param: list | None, param_name: str = "parameter if param is None: return (None, None, None, None) elif isinstance(param, list): + if len(param) == 2: # gsasii 1.0.0 RefinableParameter (A39, A84) + return (param[0], param[1], None, None) if len(param) != 4: raise ValueError( f"{param_name} must be a list of exactly 4 elements " - f"[value, refine_flag, min, max], got {len(param)} elements" + f"[value, refine_flag, min, max] or 2 elements [value, refine_flag], " + f"got {len(param)} elements" ) return tuple(param) else: @@ -224,8 +229,8 @@ def has_active_refinement_parameter(atom_param: dict) -> bool: if any(v is not None and v != [None, False, None, None] for v in value.values()): return True - elif isinstance(value, list) and len(value) == 4: - # Check RefinementParameter list [value, refine_flag, min, max] + elif isinstance(value, list) and len(value) in (2, 4): + # Check RefinementParameter list [value, refine_flag, min, max] or [value, refine_flag] if value != [None, False, None, None]: return True diff --git a/src/powderline/gateways/gsasii/project.py b/src/powderline/gateways/gsasii/project.py index 150b86b..a0aeef3 100644 --- a/src/powderline/gateways/gsasii/project.py +++ b/src/powderline/gateways/gsasii/project.py @@ -25,6 +25,7 @@ from powderline.gateways.gsasii.constraints import atom_refinement_plan, cell_refinement_plan from powderline._status import CHECK, CROSS, INFO, WARN, emoji from dataclasses import dataclass +from powderline.gateways.gsasii.schema import gsasii_space_group # Functions to add histogram with instrument parameters from recipe.xrd_data (dict with arrays) and recipe.instrument.initialization (list of dicts) @@ -331,7 +332,7 @@ def add_phase_from_cif_dict(proj: Any, cif_data: Dict[str, Any], phase_name: str pname = G2obj.MakeUniqueLabel(pname, existing_names) # Interpret space group symbol into SGData - sg = cif_data.get("space_group", "P 1") + sg = gsasii_space_group(cif_data["space_group"]) # canonical name -> GSAS-II symbol (A105, EB-41) err, SGData = G2spc.SpcGroup(sg) # Attempt normalization if initial interpretation failed @@ -365,7 +366,7 @@ def add_phase_from_cif_dict(proj: Any, cif_data: Dict[str, Any], phase_name: str # Set unit cell parameters if available unit_cell = cif_data.get("unit_cell", {}) if unit_cell: - unit_cell_constants = [unit_cell.get("a"), unit_cell.get("b"), unit_cell.get("c"), unit_cell.get("alpha"), unit_cell.get("beta"), unit_cell.get("gamma")] + unit_cell_constants = [unit_cell[k][0] for k in ("a", "b", "c", "alpha", "beta", "gamma")] # [value, refine_flag] if len(unit_cell_constants) == 6 and None not in unit_cell_constants: phase_data['General']['Cell'][1:7] = unit_cell_constants # insert a,b,c,alpha,beta,gamma try: @@ -397,10 +398,10 @@ def add_phase_from_cif_dict(proj: Any, cif_data: Dict[str, Any], phase_name: str # This code will need to be update accordingly and use "get" with defaults for missing values. E.g., don't interpret site sym or mult. for label, atom in cif_data.get("atoms", {}).items(): - x = float(atom.get("x", 0.0)) - y = float(atom.get("y", 0.0)) - z = float(atom.get("z", 0.0)) - occ = float(atom.get("occupancy", 1.0)) + x = float(atom["x"][0]) # [value, refine_flag] + y = float(atom["y"][0]) + z = float(atom["z"][0]) + occ = float(atom["occupancy"][0]) # Get ADP type for this atom (required in schema 0.22) adp_type = atom.get("ADP", None) @@ -447,16 +448,16 @@ def add_phase_from_cif_dict(proj: Any, cif_data: Dict[str, Any], phase_name: str ) # Populate U11, U22, U33, U12, U13, U23 (indices 11-16) - atom_record[11] = float(uaniso_dict["U11"]) - atom_record[12] = float(uaniso_dict["U22"]) - atom_record[13] = float(uaniso_dict["U33"]) - atom_record[14] = float(uaniso_dict["U12"]) - atom_record[15] = float(uaniso_dict["U13"]) - atom_record[16] = float(uaniso_dict["U23"]) + atom_record[11] = float(uaniso_dict["U11"][0]) + atom_record[12] = float(uaniso_dict["U22"][0]) + atom_record[13] = float(uaniso_dict["U33"][0]) + atom_record[14] = float(uaniso_dict["U12"][0]) + atom_record[15] = float(uaniso_dict["U13"][0]) + atom_record[16] = float(uaniso_dict["U23"][0]) # Set atom_record[10] to Uiso if provided, else calculate from diagonal elements if atom.get("Uiso") is not None: - atom_record[10] = float(atom["Uiso"]) + atom_record[10] = float(atom["Uiso"][0]) else: # Calculate equivalent isotropic value from diagonal elements atom_record[10] = (atom_record[11] + atom_record[12] + atom_record[13]) / 3.0 @@ -470,7 +471,7 @@ def add_phase_from_cif_dict(proj: Any, cif_data: Dict[str, Any], phase_name: str atom_record[9] = "I" uiso = atom.get("Uiso") if uiso is not None: - atom_record[10] = float(uiso) + atom_record[10] = float(uiso[0]) else: atom_record[10] = 0.0 atom_record[11:17] = [0.0]*6 # six anisotropic Uij (zeros since not used here) @@ -543,11 +544,7 @@ def add_phases_from_dict(proj: Any, hist: Any, phases_dict: dict, print_info: bo """ for phase_name, phase_info in phases_dict.items(): - structure_info = phase_info.get('structure', {}) - - # Validation check (TODO: move this upstream in phase2) - if phase_name != structure_info.get('phase_name'): - raise NameError(f"Phase name mismatch for '{phase_name}'. Structure info name is '{structure_info.get('phase_name')}'.") + structure_info = phase_info # the phase block (gsasii 1.0.0, A93); the name is the dict key (A100) # Prevent duplicate phase names - TODO: handle upstream in schema validation / recipe maker if phase_name in [p.name for p in proj.phases()]: diff --git a/src/powderline/gateways/gsasii/report_rows.py b/src/powderline/gateways/gsasii/report_rows.py new file mode 100644 index 0000000..abd08b0 --- /dev/null +++ b/src/powderline/gateways/gsasii/report_rows.py @@ -0,0 +1,34 @@ +"""GSAS-II gateway — table rows formatted exactly as ``DataFrame.to_csv`` wrote them (A83). + +The gsasii extractors build pandas DataFrames; their report files are written +by core :mod:`powderline.reports` writers (one writer for every engine). Core's +writers leave number formatting to the caller, so this formats each cell as +pandas did (``na_rep=''``, ``float_format`` on float columns only), and the +files stay byte-identical to the ones v0.1.1 wrote. +""" + +from __future__ import annotations + +import pandas as pd + + +def pandas_csv_rows(df: pd.DataFrame, float_format: str | None = None) -> list[list[str]]: + """``df``'s rows as ``to_csv(index=False, float_format=...)`` formats them.""" + columns = [] + for name in df.columns: + col = df[name] + if col.dtype.kind == "f": + fmt = (lambda v, f=float_format: f % v) if float_format else str + # numpy scalars, not Python floats: a float32 0.1 is written "0.1", as pandas does + columns.append(["" if pd.isna(v) else fmt(v) for v in col.to_numpy()]) + else: + columns.append(["" if _missing(v) else str(v) for v in col.astype(object)]) + return [list(row) for row in zip(*columns)] if columns else [] + + +def _missing(v) -> bool: + """What ``to_csv`` writes as ``na_rep``: None, NaN of any float type, ``pd.NA``, ``NaT``.""" + try: + return bool(pd.isna(v)) + except (TypeError, ValueError): # a list or array in an object column: not a scalar missing value + return False diff --git a/src/powderline/gateways/gsasii/routing.py b/src/powderline/gateways/gsasii/routing.py new file mode 100644 index 0000000..2d09deb --- /dev/null +++ b/src/powderline/gateways/gsasii/routing.py @@ -0,0 +1,37 @@ +"""GSAS-II gateway — route a validated recipe to its executor (re/04, A39). + +A native ``gsasii.rietveld`` / ``gsasii.spf`` model runs through +``executors.run_refinement``, which consumes the model directly; a 0.26.0 +``GSASII_*`` ``RecipeModel`` runs through the transient verbatim copy in +``_legacy_executor`` until re/07 removes the 0.26.0 schema. +``powderline.kicker.run_refinement`` (used by ``kicker.run``, the GSAS-II +server and the client's in-process fallback) is this function. +""" + +from __future__ import annotations + +from pathlib import Path + +from powderline.gateways.gsasii.schema import is_native_recipe, validate_recipe + + +def run_refinement(recipe, output_dir: Path, verbose: bool = False, method: str = 'server') -> dict: + """Run a validated recipe in this process; see ``executors.run_refinement`` for the result dict.""" + if is_native_recipe(recipe): + from powderline.gateways.gsasii.executors import run_refinement as run + else: + from powderline.gateways.gsasii._legacy_executor import run_refinement as run + return run(recipe, output_dir, verbose=verbose, method=method) + + +def model_from_request(recipe_dict: dict): + """The validated model for a recipe dict received by the server or the in-process fallback. + + A native recipe validates against its schema; a 0.26.0 one as before + (``RecipeModel.model_validate``, without ``validate()``'s simulation check). + Raises pydantic's ``ValidationError``. + """ + if is_native_recipe(recipe_dict): + return validate_recipe(recipe_dict) + from powderline.schema import RecipeModel + return RecipeModel.model_validate(recipe_dict) diff --git a/src/powderline/gateways/gsasii/schema.py b/src/powderline/gateways/gsasii/schema.py new file mode 100644 index 0000000..f538fde --- /dev/null +++ b/src/powderline/gateways/gsasii/schema.py @@ -0,0 +1,698 @@ +"""GSAS-II engine schemas ``gsasii.rietveld`` and ``gsasii.spf`` (engine schema 1.0.0). + +Built on the core schema (:mod:`powderline.schema_core`): a recipe is a +:class:`~powderline.schema_core.CoreRecipe`, a phase a +``Phase[RefinableParameter]`` with the gsasii fields flat in the same block +(A93). Values are in GSAS-II's native units (A4); every refinable is +``[value, refine_flag]``: gsasii 1.0.0 exposes no bounds (A84). Invalid recipes +raise pydantic's ``ValidationError`` with every problem at its location (A70). + +Engine-free: imports no GSAS-II (enforced by an import-block test), so +``validate()`` runs in the core-only environment. The gateway's runtime layer +consumes the validated model directly (A39). + +GSAS-II rules checked here (engine-behavior register entries in brackets): + +- the space group must be one GSAS-II reads as the same setting (EB-03, EB-41); +- cell refinement groups are coarser than symmetry for oblique cells, so a + refined strict subset of one is an error (A95); +- SH/L >= 0.002 (A91, EB-37); background-peak floors (A86, EB-35); Peak List + floors (A90); +- wavelength > 0 and crystallite size > 0 (GSAS-II fails otherwise; A126); +- background and Peak List peak positions lie inside the fit window (A86, A126); +- single peak fitting: the instrument profile and the peaks' own widths are + never both requested, since GSAS-II silently ignores one of them (A127, EB-53); +- simulation (``refinement_cycles == 1``) refines nothing (A88). + +Documented defaults (A119): size/strain left out = 10 um / 0, isotropic, fixed, +reported by a structured warning (A87, A107); background left out = none (one +fixed Chebyshev term 0.0, D1). +""" + +from __future__ import annotations + +from typing import Annotated, Any, Literal, Optional, Union + +import gemmi +from pydantic import Field, StrictBool, TypeAdapter, ValidationError, field_validator, model_serializer, model_validator +from pydantic_core import InitErrorDetails, PydanticCustomError + +from powderline.compat import check_engine_schema_version +from powderline.exceptions import StructuredWarning +from powderline.schema_core import ( + UNIT_ANGSTROM, + UNIT_ARBITRARY, + UNIT_DEG_2THETA, + UNIT_DEGREE, + ChebyshevBackground, + CoreFloat, + CoreInt, + CoreModel, + CoreName, + CoreRecipe, + FitRange, + Phase, + RefinableParameter, + XRDData, + absent_not_null_schema, + check_fit_range_within_data, + check_names_unique_ignoring_case, + name_keyed, +) +from powderline.symmetry import cell_tie_groups + +# --- versions (A65, A81) ---------------------------------------------------- + +ENGINE_SCHEMA_VERSION = "1.0.0" +SUPPORTED_ENGINE_SCHEMAS = "==1.0.0" +REQUIRES_CORE_SCHEMA = "==1.0.0" + +#: The gateway's schema declarations (``capabilities()`` keys, A65). +DECLARATIONS = { + "engine_schema_version": ENGINE_SCHEMA_VERSION, + "supported_engine_schemas": SUPPORTED_ENGINE_SCHEMAS, + "requires_core_schema": REQUIRES_CORE_SCHEMA, +} + +SCHEMA_NAMES = ("gsasii.rietveld", "gsasii.spf") + +# --- units (GSAS-II native, A4) ---------------------------------------------- + +UNIT_CENTIDEG = "centidegree 2theta" +UNIT_CENTIDEG2 = "centidegree^2 2theta" +UNIT_MICROMETRE = "micrometre" +UNIT_MICROSTRAIN = "microstrain" +UNIT_DIMENSIONLESS = "dimensionless" + + +def _unit(unit: str) -> dict: + return {"unit": unit} + + +def _never_null(v: Any, info) -> Any: + """A block that may be left out is never ``null`` (the A96 convention).""" + if v is None: + raise ValueError(f"{info.field_name} must not be null; leave it out instead") + return v + + +# --- space group (A105, EB-03, EB-41) ----------------------------------------- + +#: Canonical names GSAS-II 5.7.9 cannot express: ``SpcGroup`` rejects the symbol, or reads it as +#: another setting (a general position's orbit differs). Found by the orbit test over every +#: canonical name (``tests/test_gsasii_schema_engine.py``); ``:1`` origins are refused separately. +GSASII_UNREADABLE_SETTINGS = frozenset({ + "P 21212(a)", "C 2 2 21a)", "C 2 2 2a", "F 2 2 2a", "I 2 2 2a", "P 42 21 2a", "I 2 3a", + "P 21 n m", "A b a m", +}) + + +def gsasii_space_group(name: str) -> str: + """GSAS-II's symbol for a canonical (gemmi ``xhm``) space-group name. + + GSAS-II reads spaced Hermann-Mauguin symbols, hexagonal axes and origin + choice 2 by default, and rhombohedral axes from a trailing ``" R"`` + (EB-41). Raises ``ValueError`` for a setting it cannot express: origin + choice 1 (EB-03) or one of :data:`GSASII_UNREADABLE_SETTINGS`. + """ + sg = gemmi.find_spacegroup_by_name(name) + if sg is None or sg.xhm() != name: + raise ValueError(f"{name!r} is not a canonical space-group name") + if sg.ext == "1": + raise ValueError( + f"GSAS-II cannot use origin choice 1 ({name!r}); state the structure in origin choice 2 " + f"({name[:-2]}:2) (GSAS-II quirk EB-03)") + if name in GSASII_UNREADABLE_SETTINGS: + raise ValueError( + f"GSAS-II cannot express the setting {name!r} (it does not read the symbol, or reads it as " + "another setting); transform the structure to a standard setting (GSAS-II quirk EB-41)") + return sg.hm + " R" if sg.ext == "R" else sg.hm + + +def gsasii_cell_groups(space_group: str) -> tuple[tuple[str, ...], ...]: + """GSAS-II's cell refinement groups: tied or coupled parameters that refine together. + + GSAS-II varies the reciprocal metric tensor, so an oblique cell's direct + parameters cannot be held separately: monoclinic {a, c, beta} (b-unique; + {b, c, alpha} a-unique; {a, b, gamma} c-unique), and all six for triclinic + and rhombohedral axes. Otherwise they are the symmetry tie groups (A95; + ``gateways/gsasii/constraints.py`` ``cell_dof_groups``). + """ + ties = cell_tie_groups(space_group) + lengths = ("a", "b", "c") + if ties.crystal_system == "triclinic" or space_group.endswith(":R"): + return (lengths + ("alpha", "beta", "gamma"),) + if ties.crystal_system == "monoclinic": + u = ties.unique_axis + angle = {"a": "alpha", "b": "beta", "c": "gamma"}[u] + return ((u,), tuple(n for n in lengths if n != u) + (angle,)) + return tuple(g.members for g in ties.groups) + + +# --- instrument (A99, A103) -------------------------------------------------- + + +class Radiation(CoreModel): + """Radiation: GSAS-II type ``PXC`` (powder X-ray, constant wavelength) only; one wavelength.""" + + type: Literal["PXC"] = Field(description="GSAS-II instrument type; 'PXC' only (Kα doublets, TOF and " + "neutron are not supported)") + wavelength: RefinableParameter = Field(json_schema_extra=_unit(UNIT_ANGSTROM), + description="Wavelength (GSAS-II 'Lam'), > 0") + + @field_validator("wavelength") + @classmethod + def _gsasii_wavelength(cls, v: RefinableParameter) -> RefinableParameter: + if not v.value > 0: + raise ValueError(f"wavelength must be > 0, got {v.value}") + return v + + +class Geometry(CoreModel): + """Detector geometry: GSAS-II ``Bank`` and ``Azimuth`` (plain values).""" + + bank: CoreInt = Field(description="GSAS-II 'Bank'") + azimuth: CoreFloat = Field(json_schema_extra=_unit(UNIT_DEGREE), description="GSAS-II 'Azimuth'") + + +#: GSAS-II computes Bragg (Rietveld) and Peak List profiles with SH/L >= 0.002 (EB-37). +AXIAL_DIVERGENCE_MIN = 0.002 + + +class Corrections(CoreModel): + """Instrument corrections: zero shift, polarization, axial divergence (GSAS-II 'SH/L').""" + + zero_shift: RefinableParameter = Field(json_schema_extra=_unit(UNIT_CENTIDEG), description="GSAS-II 'Zero'") + polarization: RefinableParameter = Field(json_schema_extra=_unit(UNIT_DIMENSIONLESS), + description="GSAS-II 'Polariz.'") + axial_divergence: RefinableParameter = Field( + json_schema_extra=_unit(UNIT_DIMENSIONLESS), + description=f"GSAS-II 'SH/L' (asymmetry); at least {AXIAL_DIVERGENCE_MIN} (GSAS-II's own floor, A91)") + + @field_validator("axial_divergence") + @classmethod + def _sh_l_floor(cls, v: RefinableParameter) -> RefinableParameter: + if v.value < AXIAL_DIVERGENCE_MIN: + raise ValueError( + f"axial_divergence (SH/L) must be at least {AXIAL_DIVERGENCE_MIN}, got {v.value}: GSAS-II " + f"computes Bragg and Peak List profiles with SH/L >= {AXIAL_DIVERGENCE_MIN} whatever is stated " + "(GSAS-II quirk EB-37)") + return v + + +class InstrumentBroadening(CoreModel): + """Instrument profile (Caglioti U, V, W; Lorentzian X, Y, Z).""" + + U: RefinableParameter = Field(json_schema_extra=_unit(UNIT_CENTIDEG2)) + V: RefinableParameter = Field(json_schema_extra=_unit(UNIT_CENTIDEG2)) + W: RefinableParameter = Field(json_schema_extra=_unit(UNIT_CENTIDEG2)) + X: RefinableParameter = Field(json_schema_extra=_unit(UNIT_CENTIDEG)) + Y: RefinableParameter = Field(json_schema_extra=_unit(UNIT_CENTIDEG)) + Z: RefinableParameter = Field(json_schema_extra=_unit(UNIT_CENTIDEG)) + + +class Instrument(CoreModel): + """The instrument, single source (A99, A103): the gateway builds GSAS-II's instrument parameters from it.""" + + description: Optional[str] = Field(default=None, json_schema_extra=absent_not_null_schema(), + description="Free text, e.g. the beamline") + radiation: Radiation + geometry: Geometry + corrections: Corrections + broadening: InstrumentBroadening + + _not_null = field_validator("description", mode="before")(_never_null) + + @model_serializer(mode="wrap") + def _omit_absent(self, handler): + data = handler(self) + if self.description is None: + data.pop("description", None) + return data + + +def gsasii_iparm1(instrument: Instrument) -> dict: + """GSAS-II's instrument parameters (``Iparm1``) built from the instrument block (A99, A103). + + Native form ``{key: [value, value, False]}``, keys in the order of the + committed examples (N2); the setter then writes the values and refine flags. + Pure: no GSAS-II import. + """ + r, g, c, b = instrument.radiation, instrument.geometry, instrument.corrections, instrument.broadening + values = { + "Azimuth": g.azimuth, "Bank": g.bank, "Lam": r.wavelength.value, "Polariz.": c.polarization.value, + "SH/L": c.axial_divergence.value, "Type": r.type, + "U": b.U.value, "V": b.V.value, "W": b.W.value, "X": b.X.value, "Y": b.Y.value, "Z": b.Z.value, + "Zero": c.zero_shift.value, + } + return {key: [v, v, False] for key, v in values.items()} + + +# --- background (A22.6-A22.8, A86, A89) --------------------------------------- + +#: GSAS-II evaluates background peaks with these silent floors (``GSASIIpwd.py:1094-1096``, EB-35). +BACKGROUND_PEAK_FLOORS = {"intensities": 0.1, "pv_gaussian_sigma_sq": 0.01, "pv_lorentzian_gamma": 0.1} + +_PEAK_LISTS = ("positions", "intensities", "pv_gaussian_sigma_sq", "pv_lorentzian_gamma") + + +def _peak_list_problems(peaks: CoreModel, floors: dict, quirk: str) -> list[InitErrorDetails]: + """Same length for all four lists; each floored value at least its floor; errors per entry.""" + problems = [] + lengths = {name: len(getattr(peaks, name)) for name in _PEAK_LISTS} + if len(set(lengths.values())) > 1: + problems.append(InitErrorDetails( + type=PydanticCustomError("peak_lists", "all four peak lists must have the same length, got {lengths}", + {"lengths": ", ".join(f"{k}={v}" for k, v in lengths.items())}), + loc=(), input=None)) + for name, floor in floors.items(): + for i, p in enumerate(getattr(peaks, name)): + if not p.value >= floor: + problems.append(InitErrorDetails( + type=PydanticCustomError( + "below_engine_floor", + "{name}[{i}] = {value} is below {floor}, the smallest value GSAS-II computes with; " + "it would silently use {floor} (GSAS-II quirk {quirk})", + {"name": name, "i": i, "value": p.value, "floor": floor, "quirk": quirk}), + loc=(name, i), input=p.model_dump())) + return problems + + +class BackgroundPeaks(CoreModel): + """Background single peaks (pseudo-Voigt): all four lists, the same length (A22.8, A86). + + ``pv_gaussian_sigma_sq`` is the Gaussian variance (centideg², GSAS-II's + native quantity, A89); ``pv_lorentzian_gamma`` the Lorentzian width + (centideg). Each position lies inside the fit window (checked by the + payload). + """ + + positions: list[RefinableParameter] = Field(min_length=1, json_schema_extra=_unit(UNIT_DEG_2THETA)) + intensities: list[RefinableParameter] = Field(min_length=1, json_schema_extra=_unit(UNIT_ARBITRARY)) + pv_gaussian_sigma_sq: list[RefinableParameter] = Field(min_length=1, json_schema_extra=_unit(UNIT_CENTIDEG2)) + pv_lorentzian_gamma: list[RefinableParameter] = Field(min_length=1, json_schema_extra=_unit(UNIT_CENTIDEG)) + + @model_validator(mode="after") + def _lists(self) -> "BackgroundPeaks": + problems = _peak_list_problems(self, BACKGROUND_PEAK_FLOORS, "EB-35") + if problems: + raise ValidationError.from_exception_data(type(self).__name__, problems) + return self + + +def _no_background() -> ChebyshevBackground: + """Background left out = none: one fixed Chebyshev term 0.0 (D1, A119).""" + return ChebyshevBackground(num_coefficients=1, coefficients=[0.0], refine_flag=False) + + +class Background(CoreModel): + """Background: a Chebyshev polynomial (GSAS-II ``chebyschev-1``) and optional single peaks. + + Left out, ``chebyshev`` is one fixed term 0.0 (no background), replacing + GSAS-II's own constant 1.0 (EB-36, A119). Coefficients are GSAS-II's: the + polynomial is in a scaled 2theta variable and coefficient *i* multiplies + the Chebyshev polynomial of order *i* (A22.7). + """ + + chebyshev: ChebyshevBackground = Field(default_factory=_no_background) + single_peaks: Optional[BackgroundPeaks] = Field(default=None, json_schema_extra=absent_not_null_schema(), + description="Left out: no background peaks") + + _not_null = field_validator("chebyshev", "single_peaks", mode="before")(_never_null) + + @model_serializer(mode="wrap") + def _omit_absent(self, handler): + data = handler(self) + if self.single_peaks is None: + data.pop("single_peaks", None) + return data + + +# --- Peak List single peaks (gsasii.spf; A90) -------------------------------- + +#: GSAS-II evaluates Peak List peaks with these floors (``GSASIIpwd.py:1662, 1668``). +PEAK_LIST_FLOORS = {"pv_gaussian_sigma_sq": 0.001, "pv_lorentzian_gamma": 0.001} + + +class SinglePeaks(CoreModel): + """Peak List peaks for single peak fitting: all four lists, the same length. + + Widths as for background peaks (σ² in centideg², γ in centideg), each at + least GSAS-II's floor 0.001 (A90). + """ + + positions: list[RefinableParameter] = Field(min_length=1, json_schema_extra=_unit(UNIT_DEG_2THETA)) + intensities: list[RefinableParameter] = Field(min_length=1, json_schema_extra=_unit(UNIT_ARBITRARY)) + pv_gaussian_sigma_sq: list[RefinableParameter] = Field(min_length=1, json_schema_extra=_unit(UNIT_CENTIDEG2)) + pv_lorentzian_gamma: list[RefinableParameter] = Field(min_length=1, json_schema_extra=_unit(UNIT_CENTIDEG)) + + @model_validator(mode="after") + def _lists(self) -> "SinglePeaks": + problems = _peak_list_problems(self, PEAK_LIST_FLOORS, "EB-37") + if problems: + raise ValidationError.from_exception_data(type(self).__name__, problems) + return self + + +# --- phase (A93, A95, A107) --------------------------------------------------- + +#: Size/strain applied when a recipe leaves them out (A87, A107): no sample broadening. +DEFAULT_CRYSTALLITE_SIZE = 10.0 +DEFAULT_MICROSTRAIN = 0.0 +DEFAULT_LG_ETA = 1.0 + + +class SizeBroadening(CoreModel): + """Crystallite-size broadening; isotropic only in gsasii 1.0.0.""" + + model: Literal["isotropic"] = Field(description="'isotropic' (uniaxial and ellipsoidal are not supported)") + isotropic_size: RefinableParameter = Field(json_schema_extra=_unit(UNIT_MICROMETRE), description="> 0") + LG_eta: RefinableParameter = Field(json_schema_extra=_unit(UNIT_DIMENSIONLESS), + description="Lorentzian fraction of the size profile (1 = Lorentzian)") + + @field_validator("isotropic_size") + @classmethod + def _gsasii_size(cls, v: RefinableParameter) -> RefinableParameter: + if not v.value > 0: + raise ValueError(f"isotropic_size must be > 0 (a crystallite has a size), got {v.value}") + return v + + +class StrainBroadening(CoreModel): + """Microstrain broadening; isotropic only in gsasii 1.0.0.""" + + model: Literal["isotropic"] = Field(description="'isotropic' (uniaxial and generalized are not supported)") + isotropic_strain: RefinableParameter = Field(json_schema_extra=_unit(UNIT_MICROSTRAIN), + description="Microstrain, Δd/d × 10^6") + LG_eta: RefinableParameter = Field(json_schema_extra=_unit(UNIT_DIMENSIONLESS), + description="Lorentzian fraction of the strain profile (1 = Lorentzian)") + + +def _default_size() -> SizeBroadening: + return SizeBroadening(model="isotropic", isotropic_size=[DEFAULT_CRYSTALLITE_SIZE, False], + LG_eta=[DEFAULT_LG_ETA, False]) + + +def _default_strain() -> StrainBroadening: + return StrainBroadening(model="isotropic", isotropic_strain=[DEFAULT_MICROSTRAIN, False], + LG_eta=[DEFAULT_LG_ETA, False]) + + +class PeakBroadening(CoreModel): + """Sample broadening. Each part left out is 10 um / 0 microstrain, isotropic, fixed, with a warning (A107).""" + + size_broadening: SizeBroadening = Field(default_factory=_default_size) + strain_broadening: StrainBroadening = Field(default_factory=_default_strain) + + _not_null = field_validator("size_broadening", "strain_broadening", mode="before")(_never_null) + + def defaults_applied(self) -> tuple[str, ...]: + """The parts the recipe left out (the A107 default holds them).""" + return tuple(n for n in ("size_broadening", "strain_broadening") if n not in self.model_fields_set) + + +class GsasiiPhase(Phase[RefinableParameter]): + """A gsasii phase: core's space group, cell and atoms, plus ``scale`` and ``peak_broadening`` (A93). + + The phase name is the key in ``payload.phases`` (A100). GSAS-II-specific + rules: the space group must be a setting GSAS-II reads as itself (EB-03, + EB-41), and a refined cell must refine whole GSAS-II cell groups (A95). + """ + + scale: RefinableParameter = Field(json_schema_extra=_unit(UNIT_ARBITRARY), + description="Phase scale factor, >= 0 (0 allowed for reusable recipes)") + peak_broadening: PeakBroadening = Field( + default_factory=PeakBroadening, + description="Left out (in whole or part): 10 um, 0 microstrain, isotropic, fixed, with a warning (A107)") + + _gsasii_not_null = field_validator("peak_broadening", mode="before")(_never_null) + + @field_validator("scale") + @classmethod + def _gsasii_scale(cls, v: RefinableParameter) -> RefinableParameter: + if v.value < 0: + raise ValueError(f"scale must be >= 0, got {v.value}") + return v + + @field_validator("space_group") + @classmethod + def _gsasii_setting(cls, v: str) -> str: + gsasii_space_group(v) + return v + + @model_validator(mode="after") + def _gsasii_cell_groups(self) -> "GsasiiPhase": + problems = [] + for group in gsasii_cell_groups(self.space_group): + flags = [getattr(self.unit_cell, n).refine_flag for n in group] + if len(set(flags)) > 1: + refined = [n for n, f in zip(group, flags) if f] + for name, flag in zip(group, flags): + if flag != flags[0]: + problems.append(InitErrorDetails( + type=PydanticCustomError( + "gsasii_cell_group", + "GSAS-II refines {group} together (it varies the reciprocal metric tensor), but " + "only {refined} is refined; refine all of {group} or none", + {"group": ", ".join(group), "refined": ", ".join(refined)}), + loc=("unit_cell", name), input=getattr(self.unit_cell, name).model_dump())) + if problems: + raise ValidationError.from_exception_data(type(self).__name__, problems) + return self + + +# --- payloads and recipes --------------------------------------------------- + + +class RietveldControls(CoreModel): + """Refinement controls; ``refinement_cycles == 1`` is a simulation (nothing refined, A88).""" + + refinement_cycles: CoreInt = Field(default=5, ge=1, description="GSAS-II least-squares cycles") + + +class SinglePeakFittingMode(CoreModel): + use_instrument_profile: StrictBool = Field( + description="True: peak widths from the instrument profile U..Z (GSAS-II 'useIP'; the peaks' own " + "sigma^2/gamma are not used and may not be refined); False: each peak's own widths, which may " + "be refined ('hold'; U..Z have no effect and may not be refined)") + + +class SpfControls(CoreModel): + """Single-peak-fitting controls.""" + + refinement_cycles: CoreInt = Field(default=5, ge=1, description="GSAS-II peak-fit cycles") + single_peak_fitting_mode: SinglePeakFittingMode + + +def _refined_flags(obj: Any, loc: tuple = ()) -> list[tuple]: + """Locations of every refine flag set to true in a validated model tree (A88).""" + if isinstance(obj, RefinableParameter): + return [loc] if obj.refine_flag else [] + if isinstance(obj, ChebyshevBackground): + return [loc + ("refine_flag",)] if obj.refine_flag else [] + if isinstance(obj, CoreModel): + return [hit for name in type(obj).model_fields for hit in _refined_flags(getattr(obj, name), loc + (name,))] + if isinstance(obj, dict): + return [hit for k, v in obj.items() for hit in _refined_flags(v, loc + (k,))] + if isinstance(obj, list): + return [hit for i, v in enumerate(obj) for hit in _refined_flags(v, loc + (i,))] + return [] + + +class _Payload(CoreModel): + """Fields both gsasii workflows share.""" + + xrd_data: XRDData + instrument: Instrument + fit_range: Optional[FitRange] = Field(default=None, json_schema_extra=absent_not_null_schema(), + description="Left out: the data's full 2theta range") + background: Background = Field(default_factory=Background) + + _not_null = field_validator("fit_range", "background", mode="before")(_never_null) + + @model_serializer(mode="wrap") + def _omit_absent(self, handler): + data = handler(self) + if self.fit_range is None: + data.pop("fit_range", None) + return data + + @model_validator(mode="after") + def _payload_rules(self) -> "_Payload": + problems = [] + try: + check_fit_range_within_data(self.fit_range, self.xrd_data) + except ValueError as exc: + problems.append(InitErrorDetails(type=PydanticCustomError("fit_range", str(exc)), + loc=("fit_range",), input=self.fit_range.model_dump())) + if self.background.single_peaks is not None: + problems += self._outside_fit_window(self.background.single_peaks, ("background", "single_peaks"), + "background peak") + problems += self._workflow_problems() + if self.refinement_controls.refinement_cycles == 1: # simulation: nothing refined (A88) + for loc in _refined_flags(self): + problems.append(InitErrorDetails( + type=PydanticCustomError( + "simulation_refines", "refinement_cycles is 1 (a simulation), so nothing may be " + "refined; set this refine flag to false"), + loc=loc, input=True)) + if problems: + raise ValidationError.from_exception_data(type(self).__name__, problems) + return self + + def _outside_fit_window(self, peaks, loc: tuple, what: str) -> list[InitErrorDetails]: + """Peak positions outside the fit window (fit_range; open ends are the data limits; A86, A126).""" + lo, hi = self.xrd_data.tth[0], self.xrd_data.tth[-1] + if self.fit_range is not None: + lo = lo if self.fit_range.min is None else self.fit_range.min + hi = hi if self.fit_range.max is None else self.fit_range.max + return [InitErrorDetails( + type=PydanticCustomError( + "outside_fit_window", "{what} position {value} is outside the fit window " + "[{lo}, {hi}] (fit_range; open ends are the data limits)", + {"what": what, "value": p.value, "lo": lo, "hi": hi}), + loc=loc + ("positions", i), input=p.model_dump()) + for i, p in enumerate(peaks.positions) if not lo <= p.value <= hi] + + def _workflow_problems(self) -> list[InitErrorDetails]: + """Rules of one workflow's own fields (located from the payload).""" + return [] + + +class RietveldPayload(_Payload): + """``gsasii.rietveld`` payload.""" + + phases: dict[CoreName, GsasiiPhase] = Field(min_length=1, json_schema_extra=name_keyed) + refinement_controls: RietveldControls + + @field_validator("phases") + @classmethod + def _phase_names(cls, v: dict) -> dict: + check_names_unique_ignoring_case(v, "phase") + return v + + def warnings(self) -> list[StructuredWarning]: + """The size/strain defaults applied (A107), paths relative to the payload.""" + found = [] + for name, phase in self.phases.items(): + pb = phase.peak_broadening + for part in pb.defaults_applied(): + if part == "size_broadening": + what = (f"isotropic crystallite size {DEFAULT_CRYSTALLITE_SIZE:g} um", "1 um") + else: + what = (f"isotropic microstrain {DEFAULT_MICROSTRAIN:g}", "1000 microstrain") + found.append(StructuredWarning( + code="gsasii_broadening_default_applied", + message=(f"phase {name!r}: {part} left out, so PowderLine applies {what[0]}, LG_eta " + f"{DEFAULT_LG_ETA:g}, fixed (no sample broadening) instead of GSAS-II's own default " + f"({what[1]}, which broadens peaks); state {part} to choose"), + field_path=f"phases.{name}.peak_broadening.{part}", + )) + return found + + +#: Instrument profile terms that set Peak List peak widths (GSAS-II ``getCWsig`` / ``getCWgam``). +PROFILE_TERMS = ("U", "V", "W", "X", "Y", "Z") + + +class SpfPayload(_Payload): + """``gsasii.spf`` payload: single peak fitting, no phases. + + The peak widths come either from the instrument profile + (``use_instrument_profile: true``: U..Z may be refined, the peaks' own + sigma^2/gamma are not used and may not be refined) or from each peak + (``false``: the peaks' widths may be refined, the profile terms have no + effect and may not be refined). GSAS-II silently ignores a refine flag on + the side not in use (A127, EB-53), so such a flag is an error. + """ + + single_peaks: SinglePeaks + refinement_controls: SpfControls + + def _workflow_problems(self) -> list[InitErrorDetails]: + problems = self._outside_fit_window(self.single_peaks, ("single_peaks",), "Peak List peak") + use_ip = self.refinement_controls.single_peak_fitting_mode.use_instrument_profile + if use_ip: + for name in ("pv_gaussian_sigma_sq", "pv_lorentzian_gamma"): + for i, p in enumerate(getattr(self.single_peaks, name)): + if p.refine_flag: + problems.append(InitErrorDetails( + type=PydanticCustomError( + "spf_width_source", + "use_instrument_profile is true, so the peak widths come from the instrument " + "profile and GSAS-II would ignore this peak's own {name}; set its refine flag to " + "false, or set use_instrument_profile to false to refine each peak's widths " + "(GSAS-II quirk EB-53)", {"name": name}), + loc=("single_peaks", name, i), input=p.model_dump())) + else: + for name in PROFILE_TERMS: + p = getattr(self.instrument.broadening, name) + if p.refine_flag: + problems.append(InitErrorDetails( + type=PydanticCustomError( + "spf_width_source", + "use_instrument_profile is false, so each peak has its own widths and the instrument " + "profile term {name} has no effect: GSAS-II would not refine it; set its refine flag " + "to false, or set use_instrument_profile to true (GSAS-II quirk EB-53)", + {"name": name}), + loc=("instrument", "broadening", name), input=p.model_dump())) + return problems + + +class _GsasiiRecipe(CoreRecipe): + @model_validator(mode="after") + def _versions(self) -> "_GsasiiRecipe": + check_engine_schema_version(self.core_schema_version, self.engine_schema_version, + gateway="gsasii", capabilities=DECLARATIONS) + return self + + +class GsasiiRietveldRecipe(_GsasiiRecipe): + """A ``gsasii.rietveld`` recipe.""" + + schema_name: Literal["gsasii.rietveld"] + payload: RietveldPayload + + +class GsasiiSpfRecipe(_GsasiiRecipe): + """A ``gsasii.spf`` recipe.""" + + schema_name: Literal["gsasii.spf"] + payload: SpfPayload + + +GsasiiRecipe = Annotated[Union[GsasiiRietveldRecipe, GsasiiSpfRecipe], Field(discriminator="schema_name")] +_RECIPE = TypeAdapter(GsasiiRecipe, config={"title": "GsasiiRecipe"}) +_RECIPE_MODELS = {"gsasii.rietveld": GsasiiRietveldRecipe, "gsasii.spf": GsasiiSpfRecipe} + + +def is_native_recipe(recipe) -> bool: + """True for a ``gsasii.*`` recipe (dict or model); False for anything else (e.g. 0.26.0 ``GSASII_*``). + + Any ``gsasii.`` name counts, so a mistyped one (``gsasii.reitveld``) is + reported against the native schema names, not the 0.26.0 schema. + """ + if isinstance(recipe, (GsasiiRietveldRecipe, GsasiiSpfRecipe)): + return True + name = recipe.get("schema_name") if isinstance(recipe, dict) else None + return isinstance(name, str) and name.startswith("gsasii.") + + +def validate_recipe(recipe) -> GsasiiRietveldRecipe | GsasiiSpfRecipe: + """Validate a ``gsasii.*`` recipe (dict or model); return the model (idempotent, A50). + + Errors are located from the recipe's top level (``payload.phases....``); + an unknown ``schema_name`` is reported at ``schema_name``. + """ + if isinstance(recipe, (GsasiiRietveldRecipe, GsasiiSpfRecipe)): + return recipe + if not isinstance(recipe, dict): + return _RECIPE.validate_python(recipe) + model = _RECIPE_MODELS.get(recipe.get("schema_name")) + if model is None: + raise ValidationError.from_exception_data("GsasiiRecipe", [InitErrorDetails( + type=PydanticCustomError("schema_name", "schema_name '{name}' is not a gsasii workflow; expected one of " + "{names}", {"name": recipe.get("schema_name"), + "names": ", ".join(SCHEMA_NAMES)}), + loc=("schema_name",), input=recipe.get("schema_name"))]) + return model.model_validate(recipe) diff --git a/src/powderline/gateways/gsasii/server.py b/src/powderline/gateways/gsasii/server.py index e8d0816..a7a9750 100644 --- a/src/powderline/gateways/gsasii/server.py +++ b/src/powderline/gateways/gsasii/server.py @@ -87,6 +87,13 @@ class SimulationResponse(BaseModel): refined_parameters: Optional[list] = None spf_peaks: Optional[dict] = None spf_convergence_diagnostics: Optional[dict] = None + # native gsasii.* runs (re/04): standard fit statistics, GSAS-II's own values, structured warnings + r_exp: Optional[float] = None + gof: Optional[float] = None + chi2_red: Optional[float] = None + simulation_mode: Optional[bool] = None + engine_details: Optional[dict] = None + warnings: Optional[list] = None # --- Logging --- @@ -240,7 +247,7 @@ def _run_simulation(self, request: Dict[str, Any]) -> Dict[str, Any]: 'spf_convergence_diagnostics' (dict), and optional 'error', 'traceback'. """ from powderline.kicker import run_refinement - from powderline.schema import RecipeModel + from powderline.gateways.gsasii.routing import model_from_request from pydantic import ValidationError recipe_dict = request.get('recipe_data') @@ -258,7 +265,7 @@ def _run_simulation(self, request: Dict[str, Any]) -> Dict[str, Any]: # Validate recipe try: - recipe = RecipeModel.model_validate(recipe_dict) + recipe = model_from_request(recipe_dict) # native gsasii.* or 0.26.0 (re/04) except ValidationError as e: error_lines = [f"{' -> '.join(str(loc) for loc in err['loc'])}: {err['msg']}" for err in e.errors()] diff --git a/src/powderline/gateways/gsasii/setters.py b/src/powderline/gateways/gsasii/setters.py index d88f33e..1f584db 100644 --- a/src/powderline/gateways/gsasii/setters.py +++ b/src/powderline/gateways/gsasii/setters.py @@ -159,10 +159,10 @@ def set_single_peak_background(proj: Any, hist: Any, bkg_single_peaks_dict: dict proj: GSAS-II project object hist: Histogram object to set single peaks for bkg_single_peaks_dict: Dictionary with keys: - - positions: List of [[2θ, refine, min, max], ...] for peak positions - - intensities: List of [[I, refine, min, max], ...] for peak heights - - pv_gaussian_sigma: List of [[σ, refine, min, max], ...] for Gaussian widths - - pv_lorentzian_gamma: List of [[γ, refine, min, max], ...] for Lorentzian widths + - positions: List of [[2θ, refine], ...] for peak positions + - intensities: List of [[I, refine], ...] for peak heights + - pv_gaussian_sigma_sq: List of [[σ², refine], ...] Gaussian variances, centideg² (A89) + - pv_lorentzian_gamma: List of [[γ, refine], ...] for Lorentzian widths All lists must have same length (number of peaks) print_info: If True, print peak configuration to stdout @@ -175,10 +175,10 @@ def set_single_peak_background(proj: Any, hist: Any, bkg_single_peaks_dict: dict Examples: >>> # Two single peaks at 2θ=35.5° and 42.0° >>> single_peaks = { - ... 'positions': [[35.5, False, None, None], [42.0, False, None, None]], - ... 'intensities': [[50.0, True, None, None], [30.0, True, None, None]], - ... 'pv_gaussian_sigma': [[0.1, False, None, None], [0.1, False, None, None]], - ... 'pv_lorentzian_gamma': [[0.05, False, None, None], [0.05, False, None, None]] + ... 'positions': [[35.5, False], [42.0, False]], + ... 'intensities': [[50.0, True], [30.0, True]], + ... 'pv_gaussian_sigma_sq': [[0.1, False], [0.1, False]], + ... 'pv_lorentzian_gamma': [[0.15, False], [0.15, False]] ... } >>> set_single_peak_background(proj, hist, single_peaks) """ @@ -192,15 +192,12 @@ def set_single_peak_background(proj: Any, hist: Any, bkg_single_peaks_dict: dict The histogram to set the background for. bkg_single_peaks_dict : dict Dictionary containing single peak background parameters: - - positions : list of [float, bool, float, float]. Indices 2 and 3 can be None. - List of [position value, refine flag, min, max] for each peak. + - positions : list of [float, bool]: [position value, refine flag] for each peak. These are two-theta values. Q / d can be passed upstream and converted in validation. - - intensities: list of [float, bool, float, float]. Indices 2 and 3 can be None. - List of [intensity value, refine flag, min, max] for each peak. - - pv_gaussian_sigma: list of [float, bool, float, float]. Indices 2 and 3 can be None. - List of [gaussian sigma value, refine flag, min, max] for each peak. - - pv_lorentzian_gamma: list of [float, bool, float, float]. Indices 2 and 3 can be None. - List of [lorentzian gamma value, refine flag, min, max] for each peak. + - intensities: list of [float, bool]: [intensity value, refine flag] for each peak. + - pv_gaussian_sigma_sq: list of [float, bool]: [Gaussian variance (centideg²), refine flag] + for each peak. + - pv_lorentzian_gamma: list of [float, bool]: [Lorentzian gamma value, refine flag] for each peak. print_info : bool Whether to print information about the set background. Returns: @@ -209,7 +206,7 @@ def set_single_peak_background(proj: Any, hist: Any, bkg_single_peaks_dict: dict # Parse the bkg_single_peaks_dict positions = bkg_single_peaks_dict.get('positions', []) intensities = bkg_single_peaks_dict.get('intensities', []) - gaussian_sigmas = bkg_single_peaks_dict.get('pv_gaussian_sigma', []) + gaussian_sigmas = bkg_single_peaks_dict.get('pv_gaussian_sigma_sq', []) lorentzian_gammas = bkg_single_peaks_dict.get('pv_lorentzian_gamma', []) # Default intensity, sigma, gamma values to use if value input is None @@ -229,10 +226,10 @@ def set_single_peak_background(proj: Any, hist: Any, bkg_single_peaks_dict: dict single_peak_bkg = [] for i in range(num_peaks): - pos, pos_refine, pos_min, pos_max = positions[i] - inten, inten_refine, inten_min, inten_max = intensities[i] - g_sigma, g_refine, g_sigma_min, g_sigma_max = gaussian_sigmas[i] - l_gamma, l_refine, l_gamma_min, l_gamma_max = lorentzian_gammas[i] + pos, pos_refine, pos_min, pos_max = unpack_refinement_parameter(positions[i], "positions") + inten, inten_refine, inten_min, inten_max = unpack_refinement_parameter(intensities[i], "intensities") + g_sigma, g_refine, g_sigma_min, g_sigma_max = unpack_refinement_parameter(gaussian_sigmas[i], "pv_gaussian_sigma_sq") + l_gamma, l_refine, l_gamma_min, l_gamma_max = unpack_refinement_parameter(lorentzian_gammas[i], "pv_lorentzian_gamma") # TODO: move default behavior upstream # # Use default values if None provided @@ -490,7 +487,7 @@ def set_phase_unit_cell(proj: Any, phase_name: str, unit_cell_dict: dict, print_ for i, param in enumerate(cell_params): if unit_cell_dict.get(param) is not None: - value, refine_flag, min_val, max_val = unit_cell_dict[param] + value, refine_flag, min_val, max_val = unpack_refinement_parameter(unit_cell_dict[param], param) # Set the unit cell parameter value - a=1, b=2, c=3, alpha=4, beta=5, gamma=6 if value is not None: @@ -864,7 +861,7 @@ def set_phase_atom_parameters(proj: Any, phase_name: str, atom_parameters_dict: if coord in atom_params and atom_params[coord] is not None: param_value = atom_params[coord] if param_value != [None, False, None, None]: - value, refine_flag, min_val, max_val = param_value + value, refine_flag, min_val, max_val = unpack_refinement_parameter(param_value, coord) # Update coordinate value if provided if value is not None: @@ -874,7 +871,7 @@ def set_phase_atom_parameters(proj: Any, phase_name: str, atom_parameters_dict: if 'occupancy' in atom_params and atom_params['occupancy'] is not None: param_value = atom_params['occupancy'] if param_value != [None, False, None, None]: - value, refine_flag, min_val, max_val = param_value + value, refine_flag, min_val, max_val = unpack_refinement_parameter(param_value, 'occupancy') # Update occupancy value if provided if value is not None: @@ -886,7 +883,7 @@ def set_phase_atom_parameters(proj: Any, phase_name: str, atom_parameters_dict: if 'Uiso' in atom_params and atom_params['Uiso'] is not None: param_value = atom_params['Uiso'] if param_value != [None, False, None, None]: - value, refine_flag, min_val, max_val = param_value + value, refine_flag, min_val, max_val = unpack_refinement_parameter(param_value, 'Uiso') # Update Uiso value if provided if value is not None: @@ -906,7 +903,7 @@ def set_phase_atom_parameters(proj: Any, phase_name: str, atom_parameters_dict: if u_key in uaniso_dict and uaniso_dict[u_key] is not None: param_value = uaniso_dict[u_key] if param_value != [None, False, None, None]: - value, refine_flag, min_val, max_val = param_value + value, refine_flag, min_val, max_val = unpack_refinement_parameter(param_value, u_key) # Update anisotropic value if provided if value is not None: @@ -988,8 +985,8 @@ def set_phase_parameterization(proj: Any, hist: Any, phases_dict: dict, print_in # This will call on a few helper functions to set scale, unit cell, peak broadening, atom parameters, etc. # Each of these helper functions will take the proj, phase name (rather than phase object), and relevant parameterization dict - # First, get parameterization dict from phase_info (already phase-specific) - param_dict = phase_info.get('parameterization', {}) + # The phase block itself holds values and refine flags (gsasii 1.0.0, A93) + param_dict = phase_info # Keys in param dict are 'scale', 'unit_cell', 'peak_broadening', and 'atoms' for now # It is not as simple as looping over keys since some have nested dicts and others hold values or lists @@ -1062,7 +1059,7 @@ def set_instrument_parameterization(proj: Any, hist: Any, instrument_param_dict: # [0] default/starting value from iprms, [1] = value, [2] = refinement flag # Wavelength - wavelength_param = instrument_param_dict.get('wavelength', None) + wavelength_param = instrument_param_dict.get('radiation', {}).get('wavelength', None) if wavelength_param is not None: value, refine_flag, min_val, max_val = unpack_refinement_parameter(wavelength_param, "wavelength") wavelength_changed = False # track if wavelength value changed @@ -1088,7 +1085,7 @@ def set_instrument_parameterization(proj: Any, hist: Any, instrument_param_dict: # Implementing those follows the same pattern as above # Polarization - polarization_param = instrument_param_dict.get('polarization', None) + polarization_param = instrument_param_dict.get('corrections', {}).get('polarization', None) if polarization_param is not None: value, refine_flag, min_val, max_val = unpack_refinement_parameter(polarization_param, "polarization") polarization_changed = False # track if polarization value changed @@ -1148,12 +1145,14 @@ def set_instrument_parameterization(proj: Any, hist: Any, instrument_param_dict: gsas_key = 'Zero' elif param_key == 'axial_divergence': gsas_key = 'SH/L' + elif param_key == 'polarization': + continue # set above ('Polariz.'); it sits in corrections (A103) #elif param_key == 'sample_height_displacement': # irrelevant to area detector data, future work for Bragg-Brentano geometry # gsas_key = 'SampHt' else: raise ValueError(f"Unknown instrument correction parameter key '{param_key}'.") - value, refine_flag, min_val, max_val = correction_param + value, refine_flag, min_val, max_val = unpack_refinement_parameter(correction_param, param_key) correction_changed = False # Set correction parameter value diff --git a/src/powderline/gateways/gsasii/sites.py b/src/powderline/gateways/gsasii/sites.py new file mode 100644 index 0000000..103c165 --- /dev/null +++ b/src/powderline/gateways/gsasii/sites.py @@ -0,0 +1,84 @@ +"""GSAS-II gateway — check GSAS-II's reading of each atom's site after phase setup (A117). + +GSAS-II misreads some orientations of the 2-fold sites in R32, R-3m and R-3c +(``:H``): it calls them general (``'1'``) with **double** the multiplicity, so the +site weight doubles and intensities are wrong even with nothing refined (EB-40; +probe ``re04_eb40_orientation_probe.py``: 4x the intensity). Some sites it cannot +name (``'sp'``), and then refining a coordinate or ``Uaniso`` crashes its +constraint setup (EB-05). Core is right at both, so the check is a run-time error +that names an equivalent position GSAS-II reads correctly: the first member of +the atom's orbit (core's order, :func:`powderline.symmetry.orbit`) that +GSAS-II's ``SytSym`` gives core's multiplicity and a name. Never a silent +substitution: the refined coordinates would differ from the recipe's. + +Runtime layer: imports GSAS-II. +""" + +from __future__ import annotations + +from fractions import Fraction + +from GSASII import GSASIIspc as G2spc + +from powderline.exceptions import EngineExecutionError +from powderline.symmetry import orbit + + +def _coordinate(v: float) -> str: + """A coordinate as a recipe states it (A120, A121): exactly when it has a finite decimal form, else to 6 decimals.""" + frac = Fraction(v).limit_denominator(48) + if abs(float(frac) - v) < 1e-9: + d = frac.denominator + while d % 2 == 0: + d //= 2 + while d % 5 == 0: + d //= 5 + return repr(float(frac)) if d == 1 else f"{float(frac):.6f}" # 1/5 -> 0.2; 1/3 -> 0.333333 + return repr(round(v, 10)) + + +def site_problems(proj, phases) -> list[str]: + """Every atom GSAS-II reads wrongly (EB-40) or cannot name while a coordinate/Uaniso is refined (EB-05). + + ``phases`` is the validated payload's ``phases`` (name -> ``GsasiiPhase``). + """ + problems = [] + for name, phase in phases.items(): + data = proj.data["Phases"][name] + sgdata = data["General"]["SGData"] + records = {rec[0]: rec for rec in data["Atoms"]} + for label, atom in phase.atoms.items(): + gsym, gmult = records[label][7], records[label][8] + core = phase.site(label).multiplicity + refined_xyz_or_uij = any(getattr(atom, k).refine_flag for k in ("x", "y", "z")) or ( + atom.Uaniso is not None and any(p.refine_flag for _, p in atom.Uaniso)) + xyz = (atom.x.value, atom.y.value, atom.z.value) + if gmult != core: + why = (f"GSAS-II reads it as site {gsym!r} with multiplicity {gmult}, but its multiplicity is {core}: " + "the site weight would be wrong, so the intensities would be (GSAS-II quirk EB-40)") + elif gsym == "sp" and refined_xyz_or_uij: + why = ("GSAS-II cannot name its site symmetry ('sp') and fails when a coordinate or Uaniso of such a " + "site is refined (GSAS-II quirk EB-05)") + else: + continue + suggestion = next((p for p in orbit(phase.space_group, xyz) + if (lambda s: s[1] == core and s[0] != "sp")(G2spc.SytSym(list(p), sgdata)[:2])), + None) + where = f"payload.phases.{name}.atoms.{label} at ({', '.join(_coordinate(v) for v in xyz)})" + if suggestion is None: + problems.append(f"{where}: {why}; no equivalent position GSAS-II reads correctly was found") + else: + aniso = ("; its Uaniso must be transformed by the same symmetry operation (U' = R U R^T), not just " + "made site-symmetric" if atom.Uaniso is not None else "") + problems.append(f"{where}: {why}; state the equivalent position " + f"({', '.join(_coordinate(v) for v in suggestion)}) instead{aniso}") + return problems + + +def check_sites(proj, recipe) -> None: + """Raise ``EngineExecutionError`` listing every atom :func:`site_problems` finds.""" + phases = getattr(recipe.payload, "phases", None) or {} + problems = site_problems(proj, phases) + if problems: + raise EngineExecutionError("GSAS-II cannot use these atom positions as stated:\n" + + "\n".join(f" - {p}" for p in problems)) diff --git a/src/powderline/gateways/gsasii/validation.py b/src/powderline/gateways/gsasii/validation.py index dace7a4..70e9084 100644 --- a/src/powderline/gateways/gsasii/validation.py +++ b/src/powderline/gateways/gsasii/validation.py @@ -219,7 +219,8 @@ def validate(recipe: RecipeModel | dict, verbose: bool = False) -> RecipeModel: verbose: If True, print simulation-mode warnings to stdout. Returns: - Validated ``RecipeModel`` instance. + Validated ``RecipeModel`` instance; for a native ``gsasii.rietveld`` / + ``gsasii.spf`` recipe, its ``GsasiiRietveldRecipe`` / ``GsasiiSpfRecipe``. Raises: pydantic.ValidationError: If the recipe fails schema validation. @@ -238,6 +239,12 @@ def validate(recipe: RecipeModel | dict, verbose: bool = False) -> RecipeModel: """ from pydantic import ValidationError # already imported at module level, re-stated for clarity + # A native ``gsasii.*`` recipe validates against its own schema (re/04, A39); + # its simulation rule is part of that schema (A88). 0.26.0 recipes below. + from powderline.gateways.gsasii.schema import is_native_recipe, validate_recipe + if is_native_recipe(recipe): + return validate_recipe(recipe) + if not isinstance(recipe, RecipeModel): recipe = RecipeModel.model_validate(recipe) # raises ValidationError on failure diff --git a/src/powderline/kicker.py b/src/powderline/kicker.py index c723e87..91f668d 100644 --- a/src/powderline/kicker.py +++ b/src/powderline/kicker.py @@ -80,8 +80,11 @@ execute_rietveld_refinement, execute_spf_refinement, SCHEMA_EXECUTORS, - run_refinement, ) +# Native gsasii.* models -> executors.run_refinement; 0.26.0 RecipeModel -> the legacy copy (re/04) +from powderline.gateways.gsasii.routing import run_refinement # noqa: F401 +from powderline.gateways.gsasii.schema import is_native_recipe +from powderline.schema_core import collect_warnings from powderline.gateways.gsasii.extractors import ( # noqa: F401 calculate_cell_esds_from_A_matrix, extract_refined_params_from_project, @@ -295,21 +298,32 @@ def run( # 1. Validate (raises ValidationError / ValueError on failure) recipe = validate(recipe, verbose=verbose) + # A native gsasii.* recipe's structured warnings come from validating the + # user's recipe and are carried to the result (A115, A119): here for + # validate_only and the in-process run, by GSASClient on its paths + native = is_native_recipe(recipe) + warnings = collect_warnings(recipe) if native else None # 2. Validate_only short-circuit if validate_only: controls = recipe.payload.refinement_controls - return { + phases = getattr(recipe.payload, 'phases', None) + summary = { 'success': True, 'rwp': None, 'elapsed_time': 0.0, 'method': 'validate_only', 'schema_name': recipe.schema_name, - 'schema_version': recipe.schema_version, - 'phases': len(recipe.payload.phases) if recipe.payload.phases else 0, + 'phases': len(phases) if phases else 0, 'refinement_cycles': controls.refinement_cycles, 'simulation_mode': controls.refinement_cycles == 1, } + if native: + summary.update(core_schema_version=recipe.core_schema_version, + engine_schema_version=recipe.engine_schema_version, warnings=warnings) + else: + summary['schema_version'] = recipe.schema_version + return summary # 3. Ensure output directory exists output_dir = Path(output_dir) @@ -329,6 +343,13 @@ def run( auto_start_server=True, ) + if native and execution_mode == 'subprocess': # GSASClient adds them on its own paths + result['warnings'] = warnings + list(result.get('warnings') or []) + elif not native: # the server's response model carries the native-only keys as null; a 0.26.0 result never had them + for key in ('r_exp', 'gof', 'chi2_red', 'simulation_mode', 'engine_details', 'warnings'): + if key in result and result[key] is None: + del result[key] + # Normalise structured data fields to DataFrames for a consistent API. # run_refinement() and the server return JSON-serializable primitives; # converting here ensures callers always receive DataFrames regardless diff --git a/src/powderline/schema_core.py b/src/powderline/schema_core.py index 25a3331..fe2a4e0 100644 --- a/src/powderline/schema_core.py +++ b/src/powderline/schema_core.py @@ -600,7 +600,7 @@ def _absent_not_null(v: Any, info) -> Any: return v -def _absent_not_null_schema(unit: Optional[str] = None): +def absent_not_null_schema(unit: Optional[str] = None): """JSON Schema of a field that may be left out but is never null (A96). The field is ``Optional`` only so that "left out" can be held as ``None``; @@ -698,13 +698,13 @@ class Atom(CoreModel, Generic[P]): y: P = Field(description="Fractional y coordinate") z: P = Field(description="Fractional z coordinate") occupancy: P = Field(description="Site occupancy, 0 to 1 inclusive (A76)") - Multiplicity: Optional[CoreInt] = Field(default=None, json_schema_extra=_absent_not_null_schema(), + Multiplicity: Optional[CoreInt] = Field(default=None, json_schema_extra=absent_not_null_schema(), description="Site multiplicity; optional, cross-checked when stated") ADP: Literal["Uiso", "Uaniso"] = Field(description="Displacement-parameter type") - Uiso: Optional[P] = Field(default=None, json_schema_extra=_absent_not_null_schema(UNIT_ANGSTROM2), + Uiso: Optional[P] = Field(default=None, json_schema_extra=absent_not_null_schema(UNIT_ANGSTROM2), description="Isotropic ADP; required when ADP is 'Uiso', else left out") Uaniso: Optional[UanisoTensor[P]] = Field( - default=None, json_schema_extra=_absent_not_null_schema(), + default=None, json_schema_extra=absent_not_null_schema(), description="All six of U11, U22, U33, U12, U13, U23; required when ADP is 'Uaniso', else left out") _not_null = field_validator("Multiplicity", "Uiso", "Uaniso", mode="before")(_absent_not_null) @@ -1364,3 +1364,45 @@ def _engine_version(cls, v: str) -> str: @classmethod def _metadata(cls, v): return _check_metadata(v) + + +# --- warnings (A23, A115) ---------------------------------------------------- + + +def collect_warnings(model: BaseModel, path: str = "") -> list[StructuredWarning]: + """Every structured warning in a validated model tree, with full field paths. + + Each model in the tree that has a ``warnings()`` method (:meth:`Phase.warnings`, + an engine payload's documented defaults) contributes its warnings, whose + paths are relative to it; ``path`` names ``model`` itself. Called on a + validated recipe, a phase warning reads ``payload.phases..atoms.