From 88618c9c5f5c8096d6e62cf7529768deb4fc70b0 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Fri, 11 Sep 2026 00:42:35 -0400 Subject: [PATCH 01/13] docs(phys): plan incoherent CTR models --- doc/design/incoherent_ctr_models.md | 1493 +++++++++++++++++++++++++++ 1 file changed, 1493 insertions(+) create mode 100644 doc/design/incoherent_ctr_models.md diff --git a/doc/design/incoherent_ctr_models.md b/doc/design/incoherent_ctr_models.md new file mode 100644 index 0000000..8053992 --- /dev/null +++ b/doc/design/incoherent_ctr_models.md @@ -0,0 +1,1493 @@ +# Incoherent CTR models: design and implementation plan + +Status: draft, updated 2026-09-11. + +This record defines an extensible, opt-in API for incoherent CTR calculations. +The first implementation is a surface whose local height follows the existing +`PoissonProfile`, but whose laterally large height domains do not all scatter +coherently. It also defines the partial-coherence extension between the +existing coherent limit and a fully incoherent domain average. + +The design deliberately does **not** implement Sinha diffuse scattering, a +height-height power spectral density, or an off-rod diffuse line shape. It is a +domain-mixture model for kinematical CTR `F2`. Domain boundaries and scattering +that they place outside the measured CTR acceptance are outside its scope. + +## Decisions + +1. Existing `SXRDCrystal.F` and `PoissonSurface.F_uc` behavior remains fully + coherent and is the default. Merely constructing a `PoissonSurface` must not + opt a calculation into incoherent averaging. +2. An incoherent kinematical `IncoherentF2Model` consumes complete coherent + state amplitudes and returns `F2 = abs(F) ** 2`, in the squared units of the + existing kinematical structure factor. It has no `F` method because a mixed + state has no unique complex amplitude. `F2` is not detector counts and is + not reflectivity. +3. In a kinematical calculation, incoherent averaging acts on `F2` for the + **complete coherent crystal state**, `abs(F_total_n) ** 2`. It does not act + on the isolated surface correction and it does not use optical reflectivity. +4. DWBA is not supported in the first implementation. A future DWBA extension + would mix `abs(r_total_n) ** 2` for independently prepared height states; + it must not mix `abs(F_h_n) ** 2` or reuse one height-averaged optical field. +5. The public separation between coherent and incoherent models must not force + the implementation to deep-copy or reevaluate the common bulk/Film stack + for every height. `SXRDCrystal` will expose an immutable decomposition of a + kinematical evaluation, and `PoissonSurface` will expose flat-height + corrections. `PoissonHeightDomains` combines those results through the + general context API. +6. The initial partial-coherence parameter is an effective + `incoherent_fraction`, denoted by `kappa` below. It is dimensionless and lies + in `[0, 1]`. Using the coherent crystal directly is the default; + constructing an incoherent wrapper is explicit opt-in. +7. A beam coherence length is not stored on `PoissonSurface`. Surface-domain + correlation is a sample property, while beam mutual coherence, footprint, + and detector acceptance are measurement properties. A later physical + coherence kernel may combine them to calculate a point-dependent `kappa`. +8. The new module is `CTRincoherent`, not `CTRintensity`. It owns a general + `IncoherentModel` wrapper/parameter contract, the kinematical + `IncoherentF2Model` quantity contract, a registry for stable model type + names, and the Poisson height-domain implementation. +9. An `IncoherentModel` wraps one primary coherent `SXRDCrystal` and presents + the fit-parameter methods already consumed by `CTROptimizer`. + `IncoherentF2Model` adds `F2(h, k, l)`. The kinematical optimizer accepts + either the coherent crystal or an `IncoherentF2Model` as its first argument; + it must not contain Poisson-specific branches. +10. Arbitrary lists of incoherent models are not accepted because sequential + incoherent averaging has no generally valid physical meaning. A joint model + must make the correlations between its state variables explicit. +11. Height-state enumeration has a fixed, non-fit configuration setting + `exact_layer_count`, initially defaulting to 10. If the candidate + distribution contains no more than this number of layers, retain every + state. For a wider distribution, choose a contiguous retained interval + from the calculated Poisson probability masses and the profile's + cumulative tail-probability target. Do not use measured CTR values or + uncertainties to choose the forward-model support. Renormalize the + retained probabilities to sum to one and report the excluded probability + mass. `exact_layer_count` is a policy switch, not a hard cap: more than ten + states may be retained when the probability-mass criterion requires them. + Here "exact" means that every state in the profile's finite candidate + support is evaluated; it does not claim evaluation of an infinite Poisson + tail. + +## Existing coherent behavior + +`PoissonProfile.occupancy` describes cumulative material occupancy at each +structural-layer offset. `surface_occupancy` takes the difference between +successive cumulative occupancies and therefore supplies the exposed fraction +of each height state. + +`PoissonSurface.createLayers` currently assigns these fractions to coherent +domain occupancies. It also assigns the cumulative rough-Film correction and +the termination-specific replacement of exposed Film material. `F_uc` sums all +of those contributions as one complex amplitude. `SXRDCrystal.F` then adds the +bulk and every top-level component before the caller takes an absolute value. + +For flat-height amplitudes `A_n(Q)` and exposed fractions `p_n`, the current +calculation is equivalent, up to the configured Poisson tail truncation, to + +```text +A_coherent(Q) = sum_n p_n A_n(Q) +F2_coherent(Q) = abs(A_coherent(Q)) ** 2. +``` + +This is the characteristic-function form used by the classical coherent CTR +roughness treatments. Harada writes the damping factor as the squared modulus +of the Fourier sum over relative areas at each step height. Dale *et al.* use +the same squared expectation for discrete height distributions. The latter +paper calls the measured scattering incoherent in its introductory wording, +but its height-distribution factor is still the coherent amplitude average in +the terminology of this design. + +## Height-domain ensemble + +Let `p_n` be the represented probability of height state `n`, and let `A_n` be +the complete coherent amplitude of the crystal when that height state covers +the local coherent patch. Each `A_n` includes: + +- the semi-infinite bulk amplitude; +- the underlying Film and all other common components below the target + surface; +- the flat Film-height correction for state `n`; +- the termination-specific surface slab for state `n`, including its + interference with the bulk and Film; +- crystal area scaling, component weight, and supported coherent-domain + transforms. + +The limiting squared structure factors are + +```text +F2_coherent = abs(sum_n p_n A_n) ** 2 +F2_incoherent = sum_n p_n abs(A_n) ** 2. +``` + +The second expression is the large-domain limit: each coherence patch sees one +flat height and their scattering strengths add without cross terms. It is not + +```text +abs(F_bulk_and_film) ** 2 + sum_n p_n abs(F_surface_n) ** 2, +``` + +because that incorrect expression removes bulk--surface and Film--surface +interference inside each domain. + +### Partial coherence + +Use + +```text +F2_kappa = (1 - kappa) F2_coherent_live + + kappa F2_incoherent_retained, +``` + +where `kappa = 0` is fully coherent and `kappa = 1` is fully +incoherent. `incoherent_fraction` is preferred over an unqualified +`coherence` parameter so that the endpoint convention is visible at every +call site. + +For an ideal complete normalized height distribution, +`F2_coherent_live = abs(sum_n p_n A_n) ** 2` and the same interpolation can be +written + +```text +F2_kappa = F2_coherent + + kappa * ( + sum_n p_n abs(A_n) ** 2 + - abs(sum_n p_n A_n) ** 2 + ). +``` + +Equivalently in that complete-support limit, the height-state coherence +matrix is + +```text +C = (1 - kappa) outer(p, p) + kappa diag(p), +F2_kappa = A.conj() @ C @ A. +``` + +This matrix is positive semidefinite for valid probabilities and +`0 <= kappa <= 1`, so the model cannot create negative `F2`. + +This equality is exact for a complete normalized height distribution. The +finite numerical ensemble renormalizes its retained Poisson probability +masses, so the state weights used for the incoherent endpoint also form a +normalized distribution. In production, the live coherent result is the +authoritative `kappa=0` endpoint; the state-sum reconstruction is a diagnostic +which must agree within a Q-dependent truncation tolerance justified from the +excluded mass and a state-amplitude envelope, or from an extended-support +convergence check. Excluded probability alone is not an amplitude-error +bound. Because the live coherent endpoint is not replaced by the finite +state-sum reconstruction, the coherence-matrix expression above is the ideal +complete-support interpretation; the implemented partial model is the convex +interpolation of the exact live coherent endpoint and the normalized finite +incoherent endpoint. + +The interpolation follows from averaging finite coherence patches. If `f_n` +is the fraction of height `n` inside one patch, `E[f_n] = p_n`, and + +```text +Cov(f_n, f_m) = kappa * (p_n delta_nm - p_n p_m), +``` + +then `E[abs(sum_n f_n A_n) ** 2]` is exactly `F2_kappa`. For a simple patch +containing `N_eff` independent, equal-area domains, `kappa = 1 / N_eff`. +This interpretation gives the required limits without claiming that one +universal formula maps a quoted coherence length to `kappa`. + +### Coherence-length extension + +The Poisson profile is a one-point height distribution. A physical +partial-coherence calculation additionally requires the lateral two-point +correlation of the domain field. On the CTR, a reduced scalar can be written +schematically as + +```text + integral d2r K_Q(r) rho_height(r) +kappa(Q) = -----------------------------------, + integral d2r K_Q(r) +``` + +where `rho_height` is the normalized surface-domain correlation and `K_Q` +contains the incident mutual coherence, illuminated-footprint autocorrelation, +and reciprocal-space acceptance. With this convention: + +- domains much larger than the effective coherence area give `kappa -> 1`; +- many independent domains within a coherence area give `kappa -> 0`. + +The projected transverse coherence is generally anisotropic and can vary with +incidence/exit geometry. Detector acceptance and bandwidth can also make the +effective mixture reflection dependent. Therefore: + +- the first implementation accepts a scalar `incoherent_fraction`; +- a later `SurfaceCoherenceKernel` may return `kappa` for each requested point; +- domain correlation lengths belong to a sample-side domain model; +- beam coherence and acceptance belong to CTR measurement/resolution metadata; +- a CTR-only scalar fit cannot separately identify domain size and beam + coherence length. AFM, transverse scans, rocking widths, or reciprocal-space + maps are needed to constrain one of them. + +Longitudinal coherence is omitted initially. For atomic step-height +differences it will normally be much longer than the corresponding optical +path difference. If it becomes relevant, it must damp individual off-diagonal +height-state terms as a function of both `(n, m)` and scattering geometry; +that cannot in general be represented by one scalar `kappa`. + +## Why the model returns `F2`, not intensity or reflectivity + +The mathematical incoherent operation squares the complete coherent amplitude +of each state and removes the appropriate cross terms. The public quantity +must still identify which forward amplitude was squared. + +For the kinematical CTR model, `SXRDCrystal.F` is the reference-lateral-cell- +normalized structure factor. Its squared quantity is therefore + +```text +F2 = abs(F) ** 2, +``` + +in squared structure-factor units (nominally electrons squared where `F` is in +electrons). It is not observable detector counts: incident flux, footprint, +polarization/Lorentz factors, detector response, acquisition time, background, +and the fitted experimental scale are not part of `F2`. + +It is also not reflectivity. Reflectivity is a dimensionless intensity ratio +formed from the optical reflection amplitude. The quantities are: + +| Forward model | Coherent state amplitude | Domain-averaged output | +|---|---|---| +| Kinematical CTR | `F_total_n` | `F2 = sum_n p_n * abs(F_total_n) ** 2` | +| DWBA reflectivity | `r_total_n = r_0_n + r_h_n` | `R = sum_n p_n * abs(r_total_n) ** 2` | +| Isolated DWBA matrix-element diagnostic | `F_h_n` | not a reflectivity observable | + +The first implementation is specifically an incoherent `F2` model. It must +not call `SXRDCrystal.specular_reflectivity`, which is an optical profile +calculation unrelated to the kinematical `SXRDCrystal.F` path. + +`F2` must never change meaning according to an optimizer mode. A future +incoherent DWBA implementation should expose a separately named `R` boundary, +or a typed forward-result object whose quantity is explicitly `"R"`; it must +not return reflectivity from a method named `F2`. + +DWBA cannot use the proposed common-prefix optimization without further +physics. Changing the flat surface height changes the one-dimensional optical +reference profile, its Fresnel amplitude `r_0`, its internal fields, and the +atomic/reference decomposition. A physically correct DWBA domain ensemble +would prepare and evaluate each height state separately and then mix the +resulting `abs(r_total_n) ** 2`. Until that exists, combining an incoherent +surface `F2` model with DWBA must raise a clear `NotImplementedError` or +`ValueError`; silently reverting to coherent behavior is forbidden. + +## API design + +### Module and responsibility boundary + +Add `orgui.datautils.xrayutils.CTRincoherent`. Keep `SXRDCrystal` as the +coherent amplitude model and make each incoherent model a squared-structure- +factor wrapper around one coherent crystal. The module contains: + +- the abstract `IncoherentModel` wrapper/parameter contract; +- the `IncoherentF2Model` kinematical quantity contract; +- shared parameter metadata and validation; +- the model-type registry used by configuration/persistence; +- a lazy kinematical evaluation context which caches coherent components; +- reusable state-ensemble mixing helpers; and +- `PoissonHeightDomains`, the first concrete model. + +The kinematical optimizer receives either an `SXRDCrystal` or an +`IncoherentF2Model` as the existing first positional model argument. The +registry is not a service locator in the numerical hot path. It exists so a +saved type name can be converted to a class and so external packages can add +models without editing an `if/elif` chain. + +Callers that do not supply an incoherent model continue to use +`abs(crystal.F(...)) ** 2`. Thus old scripts, saved crystal files, optimizer +construction, and direct `F` calls remain coherent without a compatibility +switch. + +Add `SXRDCrystal.F2(h, k, l)` as a thin coherent convenience returning +`abs(SXRDCrystal.F(h, k, l)) ** 2`. The coherent crystal and incoherent wrapper +then share an `F2` evaluation boundary, while only the coherent crystal exposes +the phase-bearing `F` method. Existing callers may continue taking `abs(F) ** 2`. + +### Contract class + +Use an abstract base class rather than only a `Protocol`. The contract should +reuse the fit API already implemented by `SXRDCrystal` and +`CTRutil.LinearFitFunctions`, rather than inventing an `expose_parameter` +dialect. The quantity-neutral base class owns one primary coherent crystal, +model-local fit parameters, and the concatenation of both parameter blocks. +Quantity-specific subclasses define numerical evaluation methods. + +The planned public surface is: + +```python +class IncoherentModel(LinearFitFunctions, ABC): + model_type: ClassVar[str] + supported_forward_models: ClassVar[frozenset[str]] + + def __init__(self, crystal: SXRDCrystal, *, name: str = "incoherent"): ... + + @property + def coherent_model(self) -> SXRDCrystal: ... + + def addFitParameter( + self, + parameter, + limits=(-np.inf, np.inf), + **keyargs, + ) -> Parameter: ... + + @property + def fitparnames(self) -> list[str]: ... + + @property + def priors(self) -> list[object]: ... + + def getInitialParameters(self, force_recalculate=False) -> np.ndarray: ... + def getStartParamAndLimits( + self, + force_recalculate=False, + ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: ... + def setParameters(self, values: ArrayLike) -> None: ... + def getFitErrors(self) -> np.ndarray: ... + def setFitErrors(self, errors: ArrayLike | None) -> None: ... + def parameter_list(self) -> list[Parameter]: ... + def parametersToDict(self) -> dict: ... + def parametersFromDict(self, data, override_values=True) -> None: ... + def clearParameters(self) -> None: ... + def validate(self, *, forward_model: str) -> None: ... + def to_config(self) -> dict: ... + + +class IncoherentF2Model(IncoherentModel, ABC): + output_quantity: ClassVar[str] = "F2" + + def F2(self, h, k, l) -> np.ndarray: ... +``` + +The model-local parameter API follows existing usage exactly: + +```python +model.addFitParameter( + "incoherent_fraction", + limits=(0.0, 1.0), + name="rough_surface incoherent_fraction", +) +``` + +`PoissonHeightDomains` defines `parameterLookup` and `parameterLookup_inv` for +its local numerical basis, just as the existing classes do. Its initial basis +contains `incoherent_fraction`; calling `addFitParameter` promotes that basis +entry into the fitted vector. A setting therefore remains fixed unless the +caller explicitly adds it as a fit parameter. Structural/numerical policy +settings such as `surface` and `exact_layer_count` are constructor/configuration +fields, not basis entries, and never appear in the optimizer parameter vector. + +The wrapper's public parameter order is always: + +```text +[incoherent-model-local parameters] [wrapped coherent-crystal parameters] +``` + +`getStartParamAndLimits`, `getInitialParameters`, `fitparnames`, `priors`, +`parameter_list`, `setParameters`, `setFitErrors`, and `getFitErrors` all use +that same order. The base class splits and delegates the coherent tail to the +wrapped `SXRDCrystal`. This makes the wrapper directly compatible with the +optimizer's existing model contract. + +Following the existing composite-model convention, `parametersToDict` returns +both a model-local subtree and the wrapped crystal's parameter subtree; +`parametersFromDict` restores both, and `clearParameters` clears both. The +initial dictionary configuration stores only the local subtree because the +factory receives the coherent crystal separately, but that is a configuration- +layer choice rather than a different public parameter API. + +The contract has these strict invariants: + +- values, lower bounds, upper bounds, names, priors, parameter records, and + errors all have the same total length and stable order; +- local parameter semantics, default names, limits, and priors follow + `LinearFitFunctions.addFitParameter`; +- local and wrapped-crystal parameter names are unique after concatenation; +- starts lie within limits and updates validate their complete input shape + before changing either local or coherent state; +- `setFitErrors(None)` clears errors in both parameter blocks; +- after errors are cleared, `getFitErrors()` follows the existing convention + and raises `ValueError` rather than returning `None`; +- `validate` resolves referenced components by stable name on + `self.coherent_model` and checks calculation-mode capabilities; and +- `IncoherentF2Model.F2` returns a real, finite, nonnegative squared structure factor + broadcast-compatible with `h`, `k`, and `l`. It never returns `sqrt(F2)`, + detector counts, reflectivity, or a complex amplitude. + +The base class should implement the parameter methods as template methods and +delegate validation to `_validate_model(forward_model)`. +`IncoherentF2Model` implements the public `F2` validation boundary and delegates +numerical work to `_evaluate_F2(context)`. This gives third-party models one +conformance path without making them reproduce vector validation and coherent- +tail delegation. + +### Evaluation context and reusable state ensembles + +The generic model contract must not be tied to `PoissonSurface`. Its numerical +input is a lazy `KinematicIncoherentContext` containing the wrapped coherent +crystal, the requested HKL arrays, and an evaluation-local cache: + +```python +context = KinematicIncoherentContext(self.coherent_model, h, k, l) +context.coherent # lazy KinematicAmplitudeResult +context.component("rough_surface") +complete_states = context.iter_component_states( + "rough_surface", + state_evaluator, +) +``` + +`IncoherentF2Model.F2(h, k, l)` constructs this context and calls the protected +model hook. `CTROptimizer` uses the same entry point, including inside +resolution quadrature. The context centralizes reference-area scaling, +component weights, outer coherent-domain transforms, shapes, and caching; an +incoherent model must not duplicate them. + +`state_evaluator` is called by the context at every required transformed HKL; +it is not an amplitude precomputed only at the original HKL. This is necessary +because an outer coherent-domain transform changes the coordinates at which a +component state must be evaluated. The public production iterator is +`iter_component_states`: each yield contains one layer number, one normalized +probability, and the complete amplitude array for that height. One-height +streaming is the initial memory contract; batching is deferred until profiling +shows that it is needed. + +Provide an optional `CoherentStateEnsembleModel(IncoherentF2Model)` base for +the common case in which a model generates probabilities `p_n` and complete +amplitudes `A_n`. It implements validated coherent, incoherent, and coherence- +matrix sums. `PoissonHeightDomains` subclasses it. A model with different +kinematical physics may subclass `IncoherentF2Model`; a future quantity such as +reflectivity starts from the quantity-neutral `IncoherentModel` instead. + +This admits, without optimizer changes, future registered models for discrete +chemical terminations, lateral phase domains, orientational/mosaic domains, +or a physically calculated mutual-coherence matrix. Those models may expose +zero, one, or many fit parameters through the same block contract. + +The optimizer accepts one root model, not a list. For example, two rough +interfaces cannot generally be evaluated by applying two independent mixers +in sequence. Their joint states need `p_nm`, or a documented independence +factorization, and their partial coherence needs the corresponding joint +coherence matrix. Such behavior belongs in a registered joint/composite model. + +### Model registration + +Registration maps the class's single stable `model_type` persistence key to +one contract-compliant class: + +```python +@register_incoherent_model +class PoissonHeightDomains(CoherentStateEnsembleModel): + model_type = "poisson_height_domains" + + def __init__( + self, + crystal, + *, + surface, + incoherent_fraction, + exact_layer_count=10, + name="incoherent", + ): ... + ... + +model = create_incoherent_model( + "poisson_height_domains", + crystal, + surface="rough_surface", + incoherent_fraction=0.35, +) +``` + +Direct construction remains the ordinary Python API. Registration must reject +duplicate keys, abstract classes, classes not derived from `IncoherentModel`, +concrete classes without a stable `output_quantity`, and built-in key +replacement. Public discovery uses immutable `IncoherentModelInfo` records +containing the key, class, short description, supported forward models, and +`output_quantity`. `available_incoherent_models()` returns a read-only mapping +from registry key to that record, so key membership and metadata lookup are +unambiguous. +Deserialization may instantiate only registered keys; it must not import an +arbitrary class path from a file. The dictionary round-trip entry point is +`create_incoherent_model_from_config(crystal, config)`; it validates the +registered type, reconstructs the non-basis settings, and restores only the +wrapper-local parameter subtree. Full GUI/session wiring is deferred. + +Third-party registration is process-local. It need not use packaging entry +points in the first implementation. Entry-point discovery can be added later +without changing the optimizer contract. + +### `CTROptimizer` contract + +Use the existing optimizer constructor shape. The incoherent wrapper is the +forward model: + +```python +fit = CTROptimizer( + domains, + ctrs, +) +``` + +`CTROptimizer` deep-copies the complete forward-model argument once. Store that +owned object as `self.model`. For a wrapper, set `self.xtal` to +`self.model.coherent_model`; for an ordinary coherent fit, `self.model` and +`self.xtal` are the same object. This preserves the established meaning of +`optimizer.xtal`, and existing callbacks and crystal constraints continue to +receive an `SXRDCrystal`. Component selection inside the copied crystal is by +stable name, not a cached component identity. + +No `incoherent_model=` constructor keyword or setter is needed. To change +between coherent and incoherent fitting, construct the optimizer with either +the original `SXRDCrystal` or an `IncoherentF2Model` wrapping it. This keeps +one authoritative forward model and prevents mismatched crystal/model pairs. + +The complete fit-vector layout is: + +```text +[resolution] [callbacks] [subclass blocks] +[incoherent-model-local parameters] [wrapped crystal parameters] +``` + +For the base `CTROptimizer`, `subclass blocks` is empty. For +`CTROptAngleCorrection`, its existing phase/amplitude block remains in that +position. Relative ordering of every pre-existing parameter is unchanged when +the first model argument is an ordinary `SXRDCrystal`, and crystal parameters +remain the tail. The incoherent wrapper supplies its combined local/crystal +block through the methods the optimizer already calls; the optimizer does not +need Poisson-specific parameter slicing. + +At `prepareFit`, the optimizer checks each block and the concatenated vector: + +- wrapper values, bounds, names, parameter records, and priors agree in length; +- full optimizer values, bounds, and names agree in length; +- the full parameter names are unique; +- all starts satisfy bounds; +- wrapper validation succeeds for its coherent crystal; and +- the selected calculation consumes the wrapper's declared `output_quantity`; + an `F2` wrapper is rejected when DWBA is enabled. + +The optimizer's `get_parameters`, `set_parameters`, `get_bounds`, `set_errors`, +`fitparnames`, and `priors` plumbing uses `self.model`, so it sees the wrapper +as one fitted model. Physics callbacks and crystal constraints use +`self.xtal`. The wrapper performs the local/coherent parameter split +internally. For a Poisson model with one fitted fraction and two crystal +parameters, and no other optional blocks, the observable contract is: + +```text +n_parameters = 3 +fitparnames = ["rough_surface incoherent_fraction", + , ] +get_parameters() = [0.35, , ] +get_bounds()[0] = [0.0, , ] +get_bounds()[1] = [1.0, , ] +``` + +`CTROptimizer.n_parameters` is the length of this frozen vector and is valid +after `prepareFit`. `get_parameters()` is the authoritative live parameter +vector and reflects changes made to the model after preparation. At the end +of `prepareFit`, legacy `startp` is set once from `get_parameters()` so it has +the same full-vector order and length at that moment. It remains the +preparation-time snapshot and is not kept synchronized with later value +updates; callers should use `get_parameters()` instead. `startp` may be +deprecated separately in the future. Before preparation, callers can use +`len(model.getInitialParameters())` for the wrapper's total, but should not +infer the optimizer total because callbacks, subclass blocks, and resolution +settings can add prefixes. + +`prepareFit` records the model parameter-name/count/bounds signature. If +callers add, remove, rename, or re-bound fit parameters through +`optimizer.model` afterward, evaluation and parameter setters must raise that +`prepareFit()` is required; they must not slice a stale layout silently. + +`CTROptimizer.priors` retains its existing model-block scope: for an +`IncoherentF2Model`, it is the wrapper's local-plus-coherent prior list. It does +not cover resolution, callback, or angle-correction prefixes because those +APIs currently expose no priors. Extending priors to the full optimizer vector +is a separate API change and is not implied by `n_parameters`. + +Direct calculation and resolution sampling use one canonical `F2` callback. +Conversion to the stored structure-factor representation is: + +```text +predicted |F| = sqrt(resolution_operator(F2)). +``` + +Resolution is applied to `F2` before the square root. Existing helpers named +`CTRresolution.sample_intensity` and `fast_convolve_intensity` provide the +correct numerical boundary when their input is explicitly documented as +`F2`; their generic names do not change the quantity. +`sample_structure_factor` keeps its public name and `|F|` output, but should +prefer a model's `F2` method. It retains `abs(F) ** 2` as a compatibility +fallback for existing crystal-like objects which implement only `F`. + +### End-to-end usage mockup + +The following is a design mockup, not executable code yet. Assume `bulk_uc`, +`film_uc`, `termination_cells`, `h`, `k`, `l`, and `ctrs` have been constructed +by the existing APIs. + +```python +import numpy as np + +from orgui.datautils.xrayutils.CTRcalc import SXRDCrystal +from orgui.datautils.xrayutils.CTRdistributions import PoissonProfile +from orgui.datautils.xrayutils.CTRfilm import Film, PoissonSurface +from orgui.datautils.xrayutils.CTRincoherent import ( # planned module + PoissonHeightDomains, + available_incoherent_models, + create_incoherent_model_from_config, +) +from orgui.datautils.xrayutils.CTRopt import CTROptimizer + + +# This remains an ordinary coherent structural model. Constructing a +# PoissonSurface does not opt into incoherent averaging. +film = Film(film_uc, name="film") +film.basis[0] = 12.0 +surface = PoissonSurface( + termination_cells, + profile=PoissonProfile(mean_change=1.5, alpha=0.6, offset=0.0), + name="rough_surface", +) +crystal = SXRDCrystal(bulk_uc, film, surface, stacking=[0, 1]) + +F_coherent = crystal.F(h, k, l) +F2_coherent = crystal.F2(h, k, l) +np.testing.assert_allclose(F2_coherent, np.abs(F_coherent) ** 2) + + +# Explicit fully incoherent large-domain limit. Its local setting is fixed +# unless addFitParameter is called for that setting. +large_domains = PoissonHeightDomains( + crystal, + surface="rough_surface", + incoherent_fraction=1.0, + exact_layer_count=10, # default; shown here to make the policy explicit +) +assert len(large_domains.getInitialParameters()) == len(crystal.fitparnames) +F2_incoherent = large_domains.F2(h, k, l) + + +# A partially coherent model with one optimizer-visible parameter. +mixed_domains = PoissonHeightDomains( + crystal, + surface="rough_surface", + incoherent_fraction=0.35, +) +mixed_domains.addFitParameter( + "incoherent_fraction", + limits=(0.0, 1.0), + name="rough_surface incoherent_fraction", +) +assert mixed_domains.fitparnames[0] == "rough_surface incoherent_fraction" + +F2_mixed = mixed_domains.F2(h, k, l) +F_magnitude_mixed = np.sqrt(F2_mixed) # only at a |F| API boundary + + +# The wrapper is the optimizer's model. The existing model-copy operation +# copies both the wrapper and its coherent crystal. +fit = CTROptimizer( + mixed_domains, + ctrs, +) +fit.prepareFit() +assert fit.xtal is fit.model.coherent_model +assert fit.n_parameters == len(fit.get_parameters()) +prepared_start = fit.startp.copy() # legacy preparation-time snapshot +assert fit.fitparnames[-(len(crystal.fitparnames) + 1)] == ( + "rough_surface incoherent_fraction" +) +prediction = fit.flat_prediction() + + +# Registry construction is useful for a UI or dictionary configuration, but is not +# required for direct Python use. +assert "poisson_height_domains" in available_incoherent_models() +config = mixed_domains.to_config() +restored = create_incoherent_model_from_config(crystal, config) +``` + +Advanced model implementations can inspect the planned coherent decomposition +without reconstructing the crystal: + +```python +evaluation = crystal.evaluate_kinematic(h, k, l) +surface_part = next( + part for part in evaluation.components + if part.name == "rough_surface" +) +common_amplitude = evaluation.total - surface_part.amplitude + +height_states = surface.flat_domain_corrections(h, k, l) +height_states.layer_numbers +height_states.probabilities +height_states.raw_retained_probability +height_states.excluded_probability +for state in height_states.iter_states(): + state.layer_number + state.probability + state.amplitude + +# Explicit diagnostic/test materialization, not the production hot path. +all_corrections = height_states.as_array() +``` + +`PoissonHeightDomains` intentionally has no `F` method. Requesting one is a +category error because the relative phase between incoherent patches is not +defined. The initial model advertises only kinematical support, so this fails +during preparation rather than silently applying the wrong averaging: + +```python +fit.set_dwba(True) +fit.prepareFit() +# ValueError: poisson_height_domains does not support the DWBA forward model +``` + +### Persistence + +Do not put incoherent-model settings or beam coherence in `.xtal` files. +`.xtal` describes the coherent sample structure, while the effective mixture +depends on the measurement coherence and resolution. The first implementation +provides a plain dictionary contract only: + +```python +config = model.to_config() +# { +# "type": "poisson_height_domains", +# "settings": { +# "surface": "rough_surface", +# "exact_layer_count": 10, +# }, +# "parameters": , +# } + +restored = create_incoherent_model_from_config(crystal, config) +``` + +The local parameter subtree contains `basis_0` and any registered `Parameter` +records, so it is the single persisted source for `incoherent_fraction` +whether fixed or fitted. Do not duplicate the fraction under `settings`. +`exact_layer_count` is a non-fit numerical setting and therefore belongs in +`settings`. + +The factory receives an already constructed coherent crystal, resolves `type` +through the registry, wraps the crystal, and restores the stored local subtree +through the same private local helper used by the wrapper's composite +`parametersFromDict`. The public `parametersFromDict` remains a full +local-plus-coherent round trip. This increment does not choose or modify a GUI +session format; a session layer may later store this dictionary or omit it for +coherent behavior. A future physical kernel may refer separately to sample +domain-correlation parameters and measurement coherence metadata. + +## Reusable kinematical decomposition + +Strict separation of the public models does not require repeated evaluation of +the common crystal. Refactor `SXRDCrystal.F` around an immutable result: + +```python +@dataclass(frozen=True) +class KinematicComponentAmplitude: + index: int + name: str + amplitude: complex | np.ndarray + + +@dataclass(frozen=True) +class KinematicAmplitudeResult: + bulk: complex | np.ndarray + components: tuple[KinematicComponentAmplitude, ...] + total: complex | np.ndarray +``` + +Proposed public method: + +```python +result = crystal.evaluate_kinematic(h, k, l) +``` + +Its rules are: + +- input coordinates and output units are identical to `SXRDCrystal.F`; +- `bulk` includes reference-area scaling; +- every component amplitude includes its reference-area scaling, crystal + weight, and outer coherent-domain transforms; +- `total == bulk + sum(component.amplitude for component in components)`; +- returned arrays are read-only or otherwise documented as caller-owned; +- `SXRDCrystal.F` becomes a thin wrapper returning `.total`; +- the refactor must be bitwise equal where practical and numerically equal in + every existing CTR test. + +The decomposition is useful beyond this feature, but its first consumer is the +`CTRincoherent` evaluation context. For a target surface at index `s`, it +obtains the shared +amplitude once: + +```text +A_common = result.total - result.components[s].amplitude. +``` + +It then combines `A_common` with each flat-height correction of the target +surface. A raw surface correction is not yet a crystal contribution. For outer +crystal domain `d`, the state contribution must be + +```text +Delta A_n = area_scale * component_weight + * sum_d outer_domain_occupancy_d + * Delta F_n(M_d @ hkl). + +A_n = A_common + Delta A_n. +``` + +Factor the existing area/weight/domain loop into a shared evaluator used by +both `evaluate_kinematic` and the context's component-state operation. Do not +duplicate this logic in `PoissonHeightDomains`. This avoids recalculating the +bulk and ordinary Film components for every Poisson state without dropping an +established crystal-level transform. + +### Flat-height surface evaluation + +Add a `PoissonSurface` method returning coherent corrections rather than +complete crystal amplitudes: + +```python +states = surface.flat_domain_corrections( + h, + k, + l, + exact_layer_count=10, +) + +states.layer_numbers # structural-layer offsets +states.probabilities # normalized p_n for the retained Poisson support +states.raw_retained_probability +states.excluded_probability +states.excluded_lower_probability +states.excluded_upper_probability +states.iter_states() # one correction Delta F_n at a time +states.as_array() # explicit optional materialization for diagnostics +``` + +Each returned correction for state `n` contains the flat Film-height change and +the termination replacement for that height, both at occupancy one. It must +not contain `p_n`; probability is applied only by +`CoherentStateEnsembleModel`. + +For sorted height states, calculate the Film part cumulatively: + +1. Evaluate every distinct generated Film layer amplitude once. +2. Use a prefix sum from the lowest represented height to construct the flat + Film correction for every state. +3. Evaluate each termination replacement once per termination-cycle state and + translated height. +4. Add the appropriate termination replacement to each Film prefix. + +This changes the expensive scaling from repeated construction of a complete +crystal for every height to one common-crystal evaluation plus one pass over +the represented layers and height states. + +The first implementation should reuse the current layer-cycle, +growth/etching, offset, strain, termination-bank, and profile support logic. It +must not rederive layer numbering independently. It must, however, obtain the +height-state masses from `PoissonProfile.probability` rather than relying on +the final-bin convention of `surface_occupancy`; otherwise a finite-support +tail could be folded into one flat-height state before the incoherent model +applies its own cutoff policy. + +`PoissonHeightDomains.exact_layer_count` defaults to 10 and is a positive +integer, non-fit setting. If the candidate support contains at most this many +nonzero states, retain them all. For a wider support, retain a contiguous +interval containing the Poisson mode and at least `exact_layer_count` states, +then expand it using the adjacent calculated probability masses until the +profile's cumulative `tail_probability` target is met. This is a +probability-data cutoff, not a dependency on observed CTR values. It does not +impose a ten-state cap. + +Normalize the retained masses so `sum(states.probabilities) == 1` within +floating-point tolerance. Record the raw retained mass and excluded lower and +upper tail masses before normalization. Compare coherent reconstructions with +a Q-dependent tolerance derived from an amplitude envelope or an +extended-support convergence check, not directly with the dimensionless tail +probability. + +### Evaluation lifetime and cache invalidation + +The mandatory cache is an evaluation-local decomposition: for one requested +HKL batch, bulk and common components are calculated once and reused for every +height state. This removes the multiplicative height-state overhead without a +long-lived invalidation hazard. + +Do not add a global array cache in the first increment. Optimizer parameters, +component weights, coherent-domain transforms, stacking, attenuation, +reference-cell transforms, energy, and HKL arrays can all invalidate it. If +profiling later shows that cross-iteration reuse matters, add component +fingerprints or mutation generations before retaining arrays across calls. +Object identity or only `_basis_created` is not an adequate cache key. + +### Stacking boundary + +The first implementation targets one `PoissonSurface` which is the highest +material component in stacking order. This is the case in which all material +below the surface is common and only the surface correction varies. + +If a component is stacked above the target surface, its absolute position or +structure may depend on the selected height. Treating it as part of +`A_common` would be wrong. Initially reject that topology with an actionable +error. A later extension can split the evaluation into: + +- a cached common prefix below the target surface; +- the target flat-height correction; +- a state-dependent suffix restacked and evaluated for each height. + +Multiple incoherently rough interfaces also require a joint height-state +distribution or an explicit independence assumption and are out of the first +scope. + +## Validation invariants + +### Scientific endpoint tests + +1. **No opt-in, no change:** existing `SXRDCrystal.F`, optimizer predictions, + and resolution calculations remain numerically unchanged. +2. **Coherent reconstruction:** `sum(p_n * A_n)` agrees with the existing + coherent `PoissonSurface` result within a Q-dependent truncation tolerance + justified by an extended-support reference or amplitude-weighted tail. +3. **Fully incoherent endpoint:** `kappa=1` agrees with a direct weighted sum + of `F2` values from independently constructed deterministic flat-height + crystals. +4. **Deterministic height:** when only one `p_n` is nonzero, predictions are + independent of `kappa`. +5. **Two-height anti-Bragg case:** equal populations of two height states whose + amplitudes differ by a sign cancel coherently; the fully incoherent result + remains the single-flat-state `F2`, and the partial model fills the + minimum linearly with `kappa`. +6. **Nonnegative `F2`:** random valid amplitudes, probabilities, and `kappa` + values never produce negative `F2` beyond roundoff. +7. **Complete-state interference:** a fixture with a nonzero bulk amplitude + distinguishes the correct `sum(p_n * abs(A_n)**2)` result from the + incorrect sum of isolated surface squared amplitudes. + +### Existing surface behavior + +Exercise all of the current high-risk cases: + +- positive growth and negative etching; +- fractional deterministic step and nonzero `offset`; +- multi-layer Film cycles and termination-specific slabs; +- non-unit surface/Film out-of-plane repeat ratios; +- coherent-domain transforms and occupancies; +- reference-area scaling and component weights; +- Poisson probability cutoff, normalization, and tail metadata; +- copied crystals, fitted termination parameters, and coupled parameters. + +### Observable and resolution tests + +- `SXRDCrystal.F2(h, k, l)` equals `abs(SXRDCrystal.F(h, k, l)) ** 2` + and retains the reference-cell normalization of `F`. +- No resolution: stored structure-factor prediction is `sqrt(F2_kappa)`. +- Fast convolution: convolve `F2_kappa`, then take the square root. +- Quadrature sampling: evaluate `F2_kappa` at every sampled HKL, integrate, then + take the square root. +- Attempting to use the Poisson model with DWBA raises explicitly before any + mixed prediction is returned. +- Attempting to request `F` from the incoherent model is unsupported by API, + not implemented as `sqrt(F2)` with a fabricated phase. + +### Contract and optimizer tests + +- A dummy registered wrapper with zero local fit parameters contributes only + its wrapped crystal's vector, bounds, names, priors, and errors. +- Dummy models with one and several local fit parameters prove that every + contract array has the declared length and stable order. +- The complete optimizer vector has + `[resolution][callbacks][subclass][incoherent][crystal]` order, while a fit + with an ordinary coherent crystal retains its exact existing layout. +- Optimizer `set_parameters` and `set_errors` deliver the combined model block + to the wrapper; the wrapper sends only its coherent tail to `SXRDCrystal`. +- Bad shapes, duplicate names, inverted bounds, starts outside bounds, + nonfinite values, invalid errors, and missing model priors for unbounded + parameters fail during `prepareFit` with model and parameter names in the + error. +- Optimizer construction deep-copies the model. Mutating the caller's original + model afterward does not alter the optimizer-owned model. +- Registry lookup and configuration round-trips preserve settings, local + parameter order, bounds, current values, and priors supported by the existing + `Parameter.asdict` format. Duplicate or unknown type keys fail + deterministically. +- A minimal third-party test model can be registered and fitted without + changes to `CTRopt.py`, proving that the optimizer has no Poisson-specific + branch. + +## Implementation plan + +Implement this as a sequence of reviewable increments. Each increment has a +test gate; do not start optimizer integration until the coherent decomposition +and deterministic-height oracle agree. + +### Increment 0 -- characterize the existing coherent model + +Files: + +- `orgui/datautils/xrayutils/test/test_CTRcalc.py` +- optional shared fixtures in the existing xrayutils test package + +Work: + +1. Add a small Film plus `PoissonSurface` fixture with a short retained height + support and at least two distinct termination cells. +2. Construct independent deterministic flat-height crystals for every retained + height from explicit Film layers and the matching termination component. + These are the correctness oracle; they must use neither `PoissonSurface` + occupancy assembly nor the new state-extraction code. +3. Record coherent results for positive growth, negative etching, fractional + deterministic steps, nonzero offset, multi-layer cycles, component weight, + reference-area scaling, and an outer coherent-domain transform. +4. Add the equal-population, two-height anti-Bragg fixture and a fixture with a + nonzero bulk amplitude which detects accidental averaging of only the + surface correction. + +Gate: + +- Production code is unchanged. +- All current `test_CTRcalc.py` tests and the new characterization fixtures + pass. + +### Increment 1 -- add the common coherent `F2` boundary + +File: + +- `orgui/datautils/xrayutils/CTRcalc.py` + +Work: + +1. Add immutable `KinematicComponentAmplitude` and + `KinematicAmplitudeResult` records. +2. Factor the existing component area-scale, weight, and outer-domain loop into + one private evaluator. It must evaluate transformed HKL coordinates before + applying each domain occupancy. +3. Add `SXRDCrystal.evaluate_kinematic(h, k, l)`. Its result contains the + scaled bulk amplitude, one scaled amplitude per top-level component in + stable crystal order, and their total. +4. Make `SXRDCrystal.F` a thin wrapper returning the result's total. +5. Add `SXRDCrystal.F2` returning `abs(F) ** 2`, documented as squared + structure factor in the squared units of `F`, not counts or reflectivity. +6. Preserve scalar/array broadcasting, stacking, attenuation, reference frame, + and contiguous-array behavior already exercised by the CTR tests. + +Gate: + +- `evaluate_kinematic(...).total` is numerically equal to the pre-refactor + `F(...)` for all increment-0 fixtures. +- `F2(...) == abs(F(...)) ** 2` for scalar and array HKL input. +- The full `test_CTRcalc.py` suite passes before continuing. + +### Increment 2 -- factor Poisson height-state construction + +File: + +- `orgui/datautils/xrayutils/CTRfilm.py` + +Work: + +1. Extract the profile-support calculations currently embedded in + `PoissonSurface.createLayers` into one private helper used by both coherent + layer creation and flat-state evaluation. The helper owns support bounds, + cumulative material occupancy, exposed-height probabilities, Film + correction occupancy, retained masks, and tail probability metadata. Keep + the coherent `createLayers` behavior unchanged. +2. Add an immutable flat-height result containing layer numbers, + normalized probabilities, raw retained mass, excluded lower/upper mass, + `iter_states()`, and an explicit diagnostic `as_array()`. Returned arrays + must be caller-owned or read-only; the production path must not require + full amplitude materialization. +3. Add + `PoissonSurface.flat_domain_corrections(h, k, l, *, exact_layer_count=10)`. + For state `n`, the result is the occupancy-one component correction + relative to the same sharp Film boundary used by current `F_uc`; + probability is not folded into the amplitude. +4. Build each deterministic state from the material indicator for that height, + the corresponding termination slab, and the matching negative Film + termination. Reuse the current generated cells, strain, translation, + growth/etching, and termination-cycle conventions. +5. Use cumulative Film-layer amplitudes so all state corrections are produced + in one pass. Do not construct or deep-copy an `SXRDCrystal` per height. +6. Keep the calculation side-effect free with respect to persistent coherent + domain matrices. Parameter changes must be picked up through the same + synchronization path as `createLayers`. +7. Accept the wrapper's `exact_layer_count` policy value, default 10. For no + more than that many candidate states, retain all states. For wider support, + choose a contiguous, mode-containing interval from exact + `PoissonProfile.probability` masses and expand it until the profile's + cumulative tail target is met. Never use measured CTR data to choose the + support. +8. Renormalize retained probabilities exactly once, after selection. Report + raw retained and excluded tail masses, but do not use a dimensionless + probability by itself as an amplitude tolerance. Do not accidentally fold + the tail into the boundary state through `surface_occupancy` and then + renormalize it a second time. + +Gate: + +- Every returned state correction agrees with the independent deterministic + crystal oracle from increment 0 after crystal-level scaling is applied. +- The probability-weighted coherent state sum reconstructs the current + `PoissonSurface` amplitude within a Q-dependent tolerance established by an + extended-support convergence fixture or an amplitude-weighted tail bound. +- Supports below, at, and above `exact_layer_count=10` verify all-state versus + probability-cutoff selection, normalization, and reported excluded mass. +- Positive/negative growth, offsets, layer cycles, termination fits, and copied + crystals are covered. + +### Increment 3 -- implement the general incoherent wrapper contract + +Files: + +- new `orgui/datautils/xrayutils/CTRincoherent.py` +- new `orgui/datautils/xrayutils/test/test_CTRincoherent.py` + +Work: + +1. Implement quantity-neutral `IncoherentModel(LinearFitFunctions, ABC)` with + one primary `coherent_model` and a validation hook. Add + `IncoherentF2Model` with the `F2` boundary and abstract `_evaluate_F2` hook. +2. Compose the existing fit API in local-then-coherent order: + `getInitialParameters`, `getStartParamAndLimits`, `setParameters`, + `setFitErrors`, `getFitErrors`, `fitparnames`, `priors`, and + `parameter_list`. +3. Validate complete vector lengths before mutating either block. Preserve the + wrapped crystal as the tail so existing optimizer parameter assumptions + remain valid. + Treat zero-length local or coherent blocks explicitly: do not call a child + `getFitErrors()` for a block with no fit parameters, because existing + implementations may raise when no errors are set; concatenate an empty + array for that block instead. +4. Reuse `Parameter` and `addFitParameter` for model-local state. Implement + `parametersToDict`, `parametersFromDict`, and `clearParameters` as composite + operations over the local state and wrapped crystal, following Film and + other existing owning-model conventions. Keep private local-only + serialization helpers for the dictionary config factory; do not add a + parallel public parameter metadata class. +5. Implement `KinematicIncoherentContext`. It lazily evaluates the coherent + decomposition once per HKL request and applies component state evaluators + at every required outer-domain-transformed coordinate. Its production + iterator yields exactly one complete height-state amplitude at a time. +6. Add `CoherentStateEnsembleModel(IncoherentF2Model)` for validated + probability/amplitude and coherence-matrix sums. Keep the ensemble base + general; it must not import or test for `PoissonSurface`. +7. Add `register_incoherent_model`, `available_incoherent_models`, and + `create_incoherent_model`, plus + `create_incoherent_model_from_config(crystal, config)`. Use + `cls.model_type` as the only key source; reject duplicate/unknown keys, + abstract or wrong base classes, and arbitrary import paths. The config + factory restores only wrapper-local state around the supplied coherent + crystal. + +Gate: + +- Synthetic zero-, one-, and multi-local-parameter wrapper classes satisfy the + complete fit contract. +- Copying and serialization retain local parameters and the wrapped coherent + tail without identity-based component references. +- A synthetic two-state ensemble verifies coherent, incoherent, and partial + limits independently of the Poisson implementation. + +### Increment 4 -- implement `PoissonHeightDomains` + +Files: + +- `orgui/datautils/xrayutils/CTRincoherent.py` +- `orgui/datautils/xrayutils/test/test_CTRincoherent.py` + +Work: + +1. Register `PoissonHeightDomains` under `poisson_height_domains` and give it a + one-value local basis with `parameterLookup = {"incoherent_fraction": 0}`. + Initialize both `basis` and `basis_0` from the constructor value and expose + `incoherent_fraction` as a view of `basis[0]`; do not keep a second scalar + which can drift away from fitted state. + Add the positive-integer non-fit setting `exact_layer_count`, default 10. +2. Accept the coherent crystal as the first positional argument and a stable + target surface name. Validate exactly one matching `PoissonSurface`, a + finite fraction in `[0, 1]`, and the supported stacking topology. For the + first implementation, the target must be the final component in resolved + stacking order, must already be bound immediately above its Film, and may + have no same-level or later component whose position could depend on the + selected height. +3. Use the context to evaluate the common crystal amplitude once and to turn + raw flat-height corrections into complete state amplitudes with all area, + weight, and outer-domain factors applied. + Calls at outer-domain-transformed HKL must return identical layer-number and + probability vectors; only their amplitudes may depend on the transformed + coordinates. Reject inconsistent state metadata. +4. Stream one height at a time while accumulating the fully incoherent value + `F2_incoherent = sum(p_n * abs(A_n) ** 2)` without requiring an + `(N_height, N_q)` complex array. +5. Use the coherent result already cached in the context as the exact coherent + endpoint: + + ```text + F2_coherent = abs(context.coherent.total) ** 2 + F2 = (1 - kappa) * F2_coherent + + kappa * F2_incoherent. + ``` + + This is numerically identical to `coherent_model.F2(h, k, l)`, but the + implementation must use the cached total rather than call `F2` and trigger + a second bulk and Film evaluation. The independently reconstructed coherent + amplitude remains a diagnostic against the Q-dependent truncation + tolerance; it is not allowed to perturb the `kappa=0` endpoint. +6. Return finite nonnegative `float64` `F2`. Reject a request for `F` by not + defining that method. + +Gate: + +- `kappa=0` equals the wrapped coherent crystal's `F2`. +- `kappa=1` equals the deterministic-crystal `F2` average. +- Intermediate values are the convex interpolation, including the anti-Bragg + minimum. +- A one-height distribution is independent of `kappa`. +- Bulk--surface interference and outer-domain transforms match the direct + oracle. + +### Increment 5 -- integrate `F2` with resolution and `CTROptimizer` + +Files: + +- `orgui/datautils/xrayutils/CTRresolution.py` +- `orgui/datautils/xrayutils/CTRopt.py` +- `orgui/datautils/xrayutils/test/test_CTRresolution.py` +- `orgui/datautils/xrayutils/test/test_CTRopt.py` +- `orgui/datautils/xrayutils/test/test_CTRopt_dwba.py` + +Work: + +1. Change `CTRresolution.sample_structure_factor` to prefer a callable + `F2(h, k, l)` and retain `abs(F) ** 2` as the compatibility fallback. Keep + its public return value as resolution-broadened `|F|`. +2. Add one private optimizer helper for kinematical `F2`. It prefers + `model.F2(...)` and falls back to `abs(model.F(...)) ** 2` for existing + crystal-like objects and test doubles which implement only `F`. Route the + no-resolution prediction through `sqrt(helper(...))`. +3. For sampled resolution, evaluate `F2` at every quadrature point before + integration. For fast convolution, convolve central-grid `F2` directly and + take one square root afterward. +4. Continue accepting ordinary coherent crystals as the first optimizer + argument. Accept an `IncoherentF2Model` through the same argument; do not + add an `incoherent_model=` keyword or a Poisson-specific branch. +5. Store the owned forward model as `self.model`. Keep `self.xtal` as the + primary coherent `SXRDCrystal` (`self.model.coherent_model` for a wrapper), + so existing callbacks, constraints, and public `optimizer.xtal` access do + not receive a new object type. Use `self.model` for forward evaluation and + model-parameter methods. +6. Let the wrapper's existing fit methods supply local plus coherent model + parameters. Verify total values, bounds, names, priors, and error slices in + the existing optimizer order: + + ```text + [resolution] [callbacks] [subclass] + [incoherent local] [coherent crystal] + ``` + +7. Add `CTROptimizer.n_parameters`, derived from the prepared full vector. + After every parameter block has been prepared, set legacy `startp` once + from `get_parameters()`. Require `startp`, bounds, and `fitparnames` to have + exactly `n_parameters` entries at that point. Thereafter + `get_parameters()` is the authoritative current vector; `startp` remains a + preparation-time snapshot and may be deprecated in a later API change. + Record the prepared model name/count/bounds signature and reject later + structural parameter mutations until `prepareFit()` is called again. +8. Keep `CTROptimizer.priors` model-scoped, matching current behavior. Require + the wrapper's local-plus-coherent priors to match only its combined model + block and document that optional optimizer prefixes have no prior API. +9. Validate `output_quantity == "F2"` on the kinematical path. Reject an + incoherent `F2` wrapper before entering any DWBA code; never reinterpret its + result as reflectivity. +10. Preserve cache invalidation after local or coherent parameters change. + +Gate: + +- With an ordinary `SXRDCrystal`, predictions and parameter layout remain + numerically unchanged. +- Existing third-party/test crystal-like objects providing only `F` continue + to work in direct and resolution-sampled fits. +- At `prepareFit`, `startp`, bounds, `fitparnames`, `get_parameters()`, and + `n_parameters` describe the same full vector with every optional prefix + combination. Later value changes appear in `get_parameters()` without a + requirement to mutate the legacy `startp` snapshot. +- `priors` remains explicitly model-scoped and the wrapper's prior length + matches its local-plus-coherent model block. +- With a wrapper, no-resolution, sampled-resolution, and fast-convolution + predictions equal direct `F2` reference calculations. +- Resolution is demonstrably applied before `sqrt`. +- Callback, resolution, angle-correction, wrapper-local, and crystal error + slices reach the correct owners. +- Registered fit callbacks and crystal displacement constraints receive + `optimizer.xtal`, not the wrapper, while model parameter methods use + `optimizer.model`. +- DWBA plus an `F2` wrapper fails during `prepareFit` with an actionable error. + +### Increment 6 -- dictionary configuration and user documentation + +Files: + +- `orgui/datautils/xrayutils/CTRincoherent.py` +- `orgui/datautils/xrayutils/test/test_CTRincoherent.py` +- `doc/source/ctr_structure_factors.rst` +- `CHANGELOG.md` +- an example under `examples/CTR/` + +Work: + +1. Implement `to_config()` as a plain dictionary containing the registered + type, non-basis settings (`surface` and `exact_layer_count`), and the local + subtree of the wrapper's `parametersToDict` payload. The local basis is the + only saved source for `incoherent_fraction`; do not duplicate it in + settings. +2. Implement `create_incoherent_model_from_config(crystal, config)`. It wraps + the supplied coherent crystal and restores the local subtree without + replaying a second coherent parameter copy. Validate unknown keys, missing + settings, and unregistered model types. Do not add this dictionary to a GUI + or optimizer-session file format in the first implementation. +3. Document `F`, `F2`, `r`, `R`, and detector counts as distinct quantities. +4. Add the coherent/partial/incoherent Poisson example from this record and + show `addFitParameter` plus optimizer usage. +5. Add a `feat(phys)` changelog entry only when the implementation and tests + land. + +Gate: + +- Dictionary round-trip preserves wrapper type, target name, + `exact_layer_count`, fixed fraction, fitted values, limits, errors, and + priors which the existing `Parameter` serialization supports. This feature + does not broaden prior serialization. +- Old sessions and `.xtal` files remain untouched and load unchanged because + this increment does not integrate the new dictionary contract with them. + +### Deferred increment -- physical coherence kernel + +Do not include this in the first implementation. A later change may add a +sample domain-correlation model and measurement-side mutual-coherence input, +with sample-plane lengths in Angstrom and point-dependent `kappa(Q)` including +footprint and reciprocal-space acceptance. Keep scalar +`incoherent_fraction` as the explicit phenomenological model. + +### Deferred increment -- GUI/session persistence + +Choose a concrete session owner and schema only when a UI or saved optimizer +workflow needs the feature. That layer may persist the dictionary returned by +`to_config()` and interpret an absent dictionary as the existing coherent +behavior. It must continue storing the coherent structure through the +canonical `.xtal` path rather than embedding a second crystal copy. + +### Verification sequence + +Run the narrowest test after each increment, then the combined scientific +suite: + +```text +pytest orgui/datautils/xrayutils/test/test_CTRcalc.py +pytest orgui/datautils/xrayutils/test/test_CTRincoherent.py +pytest orgui/datautils/xrayutils/test/test_CTRresolution.py +pytest orgui/datautils/xrayutils/test/test_CTRopt.py +pytest orgui/datautils/xrayutils/test/test_CTRopt_dwba.py +ruff check orgui +``` + +### Definition of done + +The first implementation is complete only when: + +- coherent behavior remains the default and existing coherent tests pass; +- `F2` has one documented meaning throughout the kinematical path; +- the fully incoherent result matches independent flat-height crystals; +- the optimizer consumes the wrapper through its existing model/parameter + contract; +- resolution acts on `F2` before conversion to `|F|`; +- unsupported DWBA and stacking cases fail explicitly; +- bulk and common Film amplitudes are evaluated once per HKL batch, verified by + call-count instrumentation as well as timing; and +- no long-lived cache with incomplete invalidation is introduced. + +## Performance acceptance + +For `N_q` requested points and `N_h` represented heights, the implementation +must not evaluate the semi-infinite bulk or the common Film `N_h` times. +Expected work is approximately + +```text +one common coherent evaluation ++ one evaluation of each distinct generated layer/termination amplitude ++ O(N_h * N_q) vectorized combination and `F2` accumulation. +``` + +Benchmark the direct and quadrature-resolution paths with positive growth and +negative etching. Record peak memory as well as time. The first implementation +streams exactly one height-amplitude array at a time while accumulating +`sum(p_n * abs(A_n)**2)` and the diagnostic `sum(p_n * A_n)`; it does not +materialize an `(N_h, N_q)` complex array in the production path. Any future +batching change requires profiling evidence and must preserve the public +one-state iterator or add a separate opt-in bulk API. + +## Literature basis + +- J. Harada, [*Evaluation of the roughness of a crystal surface by X-ray + scattering. I. Theoretical considerations*](https://doi.org/10.1107/S0108767392003246), + Acta Cryst. A **48**, 764--771 (1992). The coherent CTR roughness factor is + the squared modulus of the height-population Fourier sum. +- D. Dale, A. Fleet, Y. Suzuki, and J. D. Brock, + [*X-ray scattering from real surfaces: Discrete and continuous components + of roughness*](https://doi.org/10.1103/PhysRevB.74.085419), Phys. Rev. B + **74**, 085419 (2006). Discrete height distributions enter through the + characteristic function and show the coherent anti-Bragg cancellation used + in the validation fixture. +- K. Y. C. Lee *et al.*, + [*Synchrotron X-ray study of lung surfactant-specific protein SP-B in lipid + monolayers*](https://doi.org/10.1016/S0006-3495(01)75724-4), Biophys. J. + **81**, 572--585 (2001). The appendix distinguishes coherent effective-medium + averaging from area-weighted intensity averaging when domains exceed the + beam coherence length. +- H.-J. Lee *et al.*, + [*Characterizing Pattern Structures Using X-Ray + Reflectivity*](https://www.nist.gov/publications/characterizing-pattern-structures-using-x-ray-reflectivity) + (2008). Pattern periods spanning the beam coherence scale demonstrate the + boundary of the coherent effective-medium approximation. +- T. S. Lyford, S. P. Collins, P. F. Fewster, and P. A. Thomas, + [*X-ray investigation of lateral hetero-structures of inversion domains in + LiNbO3, KTiOPO4 and KTiOAsO4*](https://doi.org/10.1107/S2053273315001503), + Acta Cryst. A **71**, 255--267 (2015). Their diffraction model combines + coherent and incoherent contributions with an empirical coefficient and + shows that effective coherence can depend on reflection, bandwidth, and + detector acceptance. +- I. A. Vartanyants and I. K. Robinson, + [*Origins of decoherence in coherent X-ray diffraction + experiments*](https://doi.org/10.1016/S0030-4018(03)01558-X), Optics + Communications **222**, 29--50 (2003). The beam is fundamentally described + by a mutual-intensity function, which need not reduce to one coherence + length. + +## Settled names + +The architectural names are decisions in this record: + +- module: `CTRincoherent`; +- wrapper/parameter contract: `IncoherentModel`; +- kinematical quantity contract: `IncoherentF2Model`; +- first implementation: `PoissonHeightDomains`; +- squared-structure-factor method: `F2(h, k, l)` on both `SXRDCrystal` and + `IncoherentF2Model`; +- optimizer construction: `CTROptimizer(incoherent_model, ctrs)`, using the + existing first positional model argument; +- registry key: `poisson_height_domains`; +- coherent decomposition: `SXRDCrystal.evaluate_kinematic`; +- flat-state evaluation: `PoissonSurface.flat_domain_corrections`; and +- evaluation context: `KinematicIncoherentContext`. + +The amplitude/`F2`/reflectivity boundary, coherent default, complete-state +averaging, DWBA exclusion, parameter-block contract, registry role, +common-prefix reuse, and ownership of coherence metadata are decisions. From 26e2bd511b816a4e255d5274e88ead5a5fc25d35 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Fri, 11 Sep 2026 13:10:57 -0400 Subject: [PATCH 02/13] docs(phys): correct incoherent CTR plan after code review Three defects would have produced wrong or breaking code: the height state index, one retention mask serving two policies, and redefining CTROptimizer.startp. --- doc/design/incoherent_ctr_models.md | 228 ++++++++++++++++++++++------ 1 file changed, 183 insertions(+), 45 deletions(-) diff --git a/doc/design/incoherent_ctr_models.md b/doc/design/incoherent_ctr_models.md index 8053992..c9b7b20 100644 --- a/doc/design/incoherent_ctr_models.md +++ b/doc/design/incoherent_ctr_models.md @@ -1,6 +1,11 @@ # Incoherent CTR models: design and implementation plan -Status: draft, updated 2026-09-11. +Status: draft, updated 2026-09-11. Revised the same day after an +implementation-readiness review against the code: the height-state index +convention, the split between coherent and incoherent retention policies, and +the treatment of `CTROptimizer.startp` were all corrected before any code was +written. Each correction is stated where it applies rather than collected in a +separate errata section. This record defines an extensible, opt-in API for incoherent CTR calculations. The first implementation is a surface whose local height follows the existing @@ -69,12 +74,30 @@ that they place outside the measured CTR acceptance are outside its scope. support is evaluated; it does not claim evaluation of an infinite Poisson tail. +12. Height states are indexed by the structural layer `n` of the top filled + layer, and the mass of that state is `PoissonProfile.probability(n + 1)`. + The `+ 1` follows from `occupancy(n) = P(H > n)`: layer `n` is the top + filled layer exactly when the signed height change equals `n + 1`. +13. `CTROptimizer.startp` keeps its existing meaning as the model-block + preparation snapshot. This feature does not redefine it into a full + optimizer vector, because that would silently change its length for every + existing fit which uses callbacks or fitted resolution. + ## Existing coherent behavior `PoissonProfile.occupancy` describes cumulative material occupancy at each -structural-layer offset. `surface_occupancy` takes the difference between -successive cumulative occupancies and therefore supplies the exposed fraction -of each height state. +structural-layer offset, `occupancy(n) = P(H > n)` for the signed height +change `H`. `surface_occupancy` takes the difference between successive +cumulative occupancies and therefore supplies the exposed fraction of each +height state. + +Two consequences fix the index convention used throughout this record. Layer +`n` is the top filled layer exactly when `H == n + 1`, so the exposed mass of +state `n` is `probability(n + 1)`, not `probability(n)`. And +`surface_occupancy` is not merely a stylistic alternative to that expression: +the two agree to floating-point roundoff on every interior bin, but +`surface_occupancy` sets its terminal bin to `occupancy[-1]` and therefore +folds the complete upper tail into the highest represented height. `PoissonSurface.createLayers` currently assigns these fractions to coherent domain occupancies. It also assigns the cumulative rough-Film correction and @@ -590,12 +613,18 @@ get_bounds()[1] = [1.0, , ] `CTROptimizer.n_parameters` is the length of this frozen vector and is valid after `prepareFit`. `get_parameters()` is the authoritative live parameter -vector and reflects changes made to the model after preparation. At the end -of `prepareFit`, legacy `startp` is set once from `get_parameters()` so it has -the same full-vector order and length at that moment. It remains the -preparation-time snapshot and is not kept synchronized with later value -updates; callers should use `get_parameters()` instead. `startp` may be -deprecated separately in the future. Before preparation, callers can use +vector and reflects changes made to the model after preparation. + +`startp` is deliberately left alone. It stays the preparation-time snapshot of +the **model block only**, exactly as `prepareFit` sets it today, and +`lower_bounds`/`higher_bounds` likewise stay the unprefixed model bounds while +`self.bounds` carries the full prefixed pair. Redefining `startp` as the full +vector would silently change its length for every existing fit which registers +a callback or fits resolution, which is a breaking change this feature has no +reason to make. For a wrapper, `startp` therefore has the wrapper's +local-plus-coherent length. Its docstring gains a note that +`get_parameters()` is the authoritative full vector; deprecating `startp` +outright is a separate API change. Before preparation, callers can use `len(model.getInitialParameters())` for the wrapper's total, but should not infer the optimizer total because callbacks, subclass blocks, and resolution settings can add prefixes. @@ -655,7 +684,10 @@ surface = PoissonSurface( profile=PoissonProfile(mean_change=1.5, alpha=0.6, offset=0.0), name="rough_surface", ) -crystal = SXRDCrystal(bulk_uc, film, surface, stacking=[0, 1]) +# Stacking levels are ordering keys and must be a numpy array; a plain list +# fails when the constructor reorders them. The values match the convention +# used by the existing tests for a Film plus surface pair. +crystal = SXRDCrystal(bulk_uc, film, surface, stacking=np.array([1, 2])) F_coherent = crystal.F(h, k, l) F2_coherent = crystal.F2(h, k, l) @@ -895,10 +927,24 @@ the represented layers and height states. The first implementation should reuse the current layer-cycle, growth/etching, offset, strain, termination-bank, and profile support logic. It must not rederive layer numbering independently. It must, however, obtain the -height-state masses from `PoissonProfile.probability` rather than relying on -the final-bin convention of `surface_occupancy`; otherwise a finite-support -tail could be folded into one flat-height state before the incoherent model -applies its own cutoff policy. +height-state masses as `PoissonProfile.probability(layer_numbers + 1)` rather +than from `surface_occupancy`; otherwise the finite-support tail is folded +into the highest flat-height state before the incoherent model applies its own +cutoff policy. At `kappa = 1` that folded bin would become a single spurious +flat state carrying the whole tail mass, which is a physically wrong result +that still looks plausible. Because the interior bins of the two expressions +agree to roundoff, a regression test must assert both halves of the relation: +interior agreement, and deliberate disagreement at the terminal bin. + +The retained coherent component set and the retained height-state ensemble are +different objects and must not be conflated. `createLayers` retains the union +of exposed-surface and Film-correction layers and never renormalizes; the +incoherent ensemble retains exposed-height states only, and renormalizes. A +third range is independent of both: the Film correction of flat state `n` +covers every structural layer between the sharp boundary and `n`, so the +cumulative Film prefix must span +`[min(0, lowest_retained_state), highest_retained_state]` taken from the full +candidate support rather than from either retained mask. `PoissonHeightDomains.exact_layer_count` defaults to 10 and is a positive integer, non-fit setting. If the candidate support contains at most this many @@ -1030,12 +1076,49 @@ Implement this as a sequence of reviewable increments. Each increment has a test gate; do not start optimizer integration until the coherent decomposition and deterministic-height oracle agree. +### Preparatory fix -- error routing in `evaluateStatistics` + +Files: + +- `orgui/datautils/xrayutils/CTRopt.py` +- `orgui/datautils/xrayutils/test/test_CTRopt.py` + +Work: + +1. `CTROptimizer.evaluateStatistics` calls `self.xtal.setFitErrors(errors[3:])` + or `self.xtal.setFitErrors(errors)` directly. That hardcodes the resolution + prefix width and skips registered callbacks entirely, so with a callback + present the crystal receives the callback's error slice. Replace both + branches with `self.set_errors(errors)`, which is the existing splitter for + resolution, callbacks, subclass blocks, and the model. +2. This is a pre-existing defect independent of incoherent models, so it lands + as its own commit ahead of increment 1: bisectable on its own, and the + wrapper work does not carry an unrelated behavior change. +3. It also composes with increment 5 at no extra cost. Once `_set_model_errors` + forwards to `self.model`, a wrapper's combined local-plus-coherent error + block reaches the wrapper through the same splitter, with no second site to + update. +4. The repository AGENTS.md places error propagation under the `phys` scope, so + the commit is `fix(phys)`. The method is already `DeprecationWarning`-marked + and stays that way. + +Gate: + +- A fit with a registered callback and fitted resolution routes each error + slice to its owner, asserted per block. +- With no callbacks and no fitted resolution, `evaluateStatistics` behavior is + numerically unchanged. + ### Increment 0 -- characterize the existing coherent model Files: - `orgui/datautils/xrayutils/test/test_CTRcalc.py` -- optional shared fixtures in the existing xrayutils test package +- new `orgui/datautils/xrayutils/test/_poisson_oracle.py` for the shared + deterministic flat-height fixtures. The repository has no `conftest.py` and + its CTR tests are `unittest`-style classes run under pytest, so the shared + fixtures are a plain importable module in the test package, not pytest + fixtures. Work: @@ -1089,18 +1172,30 @@ Gate: ### Increment 2 -- factor Poisson height-state construction -File: +Files: - `orgui/datautils/xrayutils/CTRfilm.py` +- new `orgui/datautils/xrayutils/test/test_CTRfilm.py`, importing the + deterministic oracle from `test/_poisson_oracle.py` Work: 1. Extract the profile-support calculations currently embedded in - `PoissonSurface.createLayers` into one private helper used by both coherent - layer creation and flat-state evaluation. The helper owns support bounds, - cumulative material occupancy, exposed-height probabilities, Film - correction occupancy, retained masks, and tail probability metadata. Keep - the coherent `createLayers` behavior unchanged. + `PoissonSurface.createLayers` into one private helper which applies **no** + retention mask. It returns a frozen candidate record: support bounds, + layer numbers, cumulative material occupancy, exposed-height probabilities + (`probability(layer_numbers + 1)`), Film correction occupancy, and the tail + probability. Selection is then two separate callers of that record: + + - the coherent selector reproduces today's union mask + `(exposed > tail_probability) | (abs(film_correction) > tail_probability)` + with no renormalization, so `createLayers` behavior is unchanged; and + - the incoherent selector applies the `exact_layer_count` policy below to + the exposed-height states alone and renormalizes. + + Do not give the helper a single "retained mask" output. The two policies + retain different sets for different reasons, and merging them is what would + silently make the coherent path adopt the incoherent cutoff or the reverse. 2. Add an immutable flat-height result containing layer numbers, normalized probabilities, raw retained mass, excluded lower/upper mass, `iter_states()`, and an explicit diagnostic `as_array()`. Returned arrays @@ -1116,7 +1211,11 @@ Work: termination. Reuse the current generated cells, strain, translation, growth/etching, and termination-cycle conventions. 5. Use cumulative Film-layer amplitudes so all state corrections are produced - in one pass. Do not construct or deep-copy an `SXRDCrystal` per height. + in one pass. Do not construct or deep-copy an `SXRDCrystal` per height. The + prefix range is taken from the candidate support, not from either retained + mask: it must cover every structural layer between the sharp Film boundary + and the highest retained state, including layers whose own exposed + probability is negligible. 6. Keep the calculation side-effect free with respect to persistent coherent domain matrices. Parameter changes must be picked up through the same synchronization path as `createLayers`. @@ -1124,8 +1223,11 @@ Work: more than that many candidate states, retain all states. For wider support, choose a contiguous, mode-containing interval from exact `PoissonProfile.probability` masses and expand it until the profile's - cumulative tail target is met. Never use measured CTR data to choose the - support. + cumulative tail target is met. The mode is `argmax` of that calculated + probability array with ties resolved to the lower index, never + `floor(rate)`: for `alpha < 1` the distribution is a two-point step mixture + convolved with the Poisson and its maximum need not sit at the Poisson + mode. Never use measured CTR data to choose the support. 8. Renormalize retained probabilities exactly once, after selection. Report raw retained and excluded tail masses, but do not use a dimensionless probability by itself as an amplitude tolerance. Do not accidentally fold @@ -1136,6 +1238,11 @@ Gate: - Every returned state correction agrees with the independent deterministic crystal oracle from increment 0 after crystal-level scaling is applied. +- Interior bins of `surface_occupancy(n)` equal `probability(n + 1)` to + roundoff, and the terminal bin deliberately does not. Both halves are + asserted so a later refactor cannot silently exchange the two expressions. +- The coherent and incoherent selectors retain their own sets from one shared + candidate record, and `createLayers` output is unchanged by the refactor. - The probability-weighted coherent state sum reconstructs the current `PoissonSurface` amplitude within a Q-dependent tolerance established by an extended-support convergence fixture or an amplitude-weighted tail bound. @@ -1156,10 +1263,25 @@ Work: 1. Implement quantity-neutral `IncoherentModel(LinearFitFunctions, ABC)` with one primary `coherent_model` and a validation hook. Add `IncoherentF2Model` with the `F2` boundary and abstract `_evaluate_F2` hook. -2. Compose the existing fit API in local-then-coherent order: - `getInitialParameters`, `getStartParamAndLimits`, `setParameters`, - `setFitErrors`, `getFitErrors`, `fitparnames`, `priors`, and - `parameter_list`. +2. Compose the existing fit API in local-then-coherent order. Call + `super().__init__()` so the inherited empty `basis`, `basis_0` and + `parameters` exist even for a wrapper with no local parameters. Note that + `LinearFitFunctions` already implements `getInitialParameters`, + `getStartParamAndLimits`, `setFitErrors`, `getFitErrors`, `fitparnames`, + `priors` and `parameter_list` as **local-only** methods. They must + therefore all be overridden as composites which take the local block from + an explicit `LinearFitFunctions.(self, ...)` call and append the + coherent tail. Inheriting one of them by accident does not raise; it + returns a silently short vector, which is the failure this contract exists + to prevent. Only `setParameters` and `validate` are genuinely new: + `setParameters` exists on `SXRDCrystal`, which is not a + `LinearFitFunctions` subclass, and has no base implementation to reuse. + Override the inherited `setFitParameters` to raise `NotImplementedError` + naming `setParameters` as the composite entry point. Nothing calls it on a + wrapper -- `SXRDCrystal.setParameters` only calls `setFitParameters` on its + own component unit cells, and a wrapper is never a component -- so raising + costs nothing and converts a silent wrong-length write into an immediate + error. 3. Validate complete vector lengths before mutating either block. Preserve the wrapped crystal as the tail so existing optimizer parameter assumptions remain valid. @@ -1191,7 +1313,12 @@ Work: Gate: - Synthetic zero-, one-, and multi-local-parameter wrapper classes satisfy the - complete fit contract. + complete fit contract. One conformance test catches every forgotten + composite override at once: `fitparnames`, `getInitialParameters()`, + `parameter_list()`, `priors` and `getStartParamAndLimits()[0]` all have the + same length, and the tail of each equals the wrapped crystal's own entries + in order. +- `setFitParameters` raises on a wrapper instead of writing the local block. - Copying and serialization retain local parameters and the wrapped coherent tail without identity-based component references. - A synthetic two-state ensemble verifies coherent, incoherent, and partial @@ -1209,8 +1336,13 @@ Work: 1. Register `PoissonHeightDomains` under `poisson_height_domains` and give it a one-value local basis with `parameterLookup = {"incoherent_fraction": 0}`. Initialize both `basis` and `basis_0` from the constructor value and expose - `incoherent_fraction` as a view of `basis[0]`; do not keep a second scalar - which can drift away from fitted state. + `incoherent_fraction` as a property whose getter returns + `float(self.basis[0])` and whose setter validates a finite value in + `[0, 1]` before writing `self.basis[0]`. A numpy scalar view is not + available for this; the point is that there is no second stored scalar + which can drift away from fitted state. Seeding `basis_0` as well is what + makes `updateFromParameters` restore the configured fraction when it is not + a fit parameter. Add the positive-integer non-fit setting `exact_layer_count`, default 10. 2. Accept the coherent crystal as the first positional argument and a stable target surface name. Validate exactly one matching `PoissonSurface`, a @@ -1275,8 +1407,13 @@ Work: crystal-like objects and test doubles which implement only `F`. Route the no-resolution prediction through `sqrt(helper(...))`. 3. For sampled resolution, evaluate `F2` at every quadrature point before - integration. For fast convolution, convolve central-grid `F2` directly and - take one square root afterward. + integration. For fast convolution, fill the cached input collection with + `sqrt(F2)` and keep calling `CTRresolution.fast_convolve`, which squares + its input, convolves, and takes one square root. That is exactly + "resolution applied to `F2` before the square root" and it introduces no + second collection-building path, so `_require_structure_factors` and + `preserve_measurement_metadata` keep working unchanged. The cost is one + redundant square-root/square round trip per point, at roundoff level. 4. Continue accepting ordinary coherent crystals as the first optimizer argument. Accept an `IncoherentF2Model` through the same argument; do not add an `incoherent_model=` keyword or a Poisson-specific branch. @@ -1294,14 +1431,13 @@ Work: [incoherent local] [coherent crystal] ``` -7. Add `CTROptimizer.n_parameters`, derived from the prepared full vector. - After every parameter block has been prepared, set legacy `startp` once - from `get_parameters()`. Require `startp`, bounds, and `fitparnames` to have - exactly `n_parameters` entries at that point. Thereafter - `get_parameters()` is the authoritative current vector; `startp` remains a - preparation-time snapshot and may be deprecated in a later API change. - Record the prepared model name/count/bounds signature and reject later - structural parameter mutations until `prepareFit()` is called again. +7. Add `CTROptimizer.n_parameters`, derived from the prepared full vector, so + that `n_parameters == len(get_parameters()) == len(fitparnames) == + len(bounds[0])`. Leave `startp`, `lower_bounds` and `higher_bounds` with + their existing model-block scope and length; only their docstrings change, + to point at `get_parameters()` as the authoritative full vector. Record the + prepared model name/count/bounds signature and reject later structural + parameter mutations until `prepareFit()` is called again. 8. Keep `CTROptimizer.priors` model-scoped, matching current behavior. Require the wrapper's local-plus-coherent priors to match only its combined model block and document that optional optimizer prefixes have no prior API. @@ -1316,10 +1452,11 @@ Gate: numerically unchanged. - Existing third-party/test crystal-like objects providing only `F` continue to work in direct and resolution-sampled fits. -- At `prepareFit`, `startp`, bounds, `fitparnames`, `get_parameters()`, and +- At `prepareFit`, `bounds`, `fitparnames`, `get_parameters()` and `n_parameters` describe the same full vector with every optional prefix - combination. Later value changes appear in `get_parameters()` without a - requirement to mutate the legacy `startp` snapshot. + combination. `startp` is explicitly excluded from that equality and instead + satisfies `len(startp) == len(model.getInitialParameters())`, unchanged for + an ordinary coherent fit with callbacks or fitted resolution. - `priors` remains explicitly model-scoped and the wrapper's prior length matches its local-plus-coherent model block. - With a wrapper, no-resolution, sampled-resolution, and fast-convolution @@ -1392,6 +1529,7 @@ suite: ```text pytest orgui/datautils/xrayutils/test/test_CTRcalc.py +pytest orgui/datautils/xrayutils/test/test_CTRfilm.py pytest orgui/datautils/xrayutils/test/test_CTRincoherent.py pytest orgui/datautils/xrayutils/test/test_CTRresolution.py pytest orgui/datautils/xrayutils/test/test_CTRopt.py From a2ff4510e54fc4365e45da1fe2fd07b4f669cff6 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Fri, 11 Sep 2026 13:16:48 -0400 Subject: [PATCH 03/13] docs(phys): resolve remaining inconsistencies in the CTR plan The symbol p_n meant both the coherent exposed fractions and the normalized ensemble masses; the coherent ones are now q_n. Also fixes the Film prefix range and the mode definition. --- doc/design/incoherent_ctr_models.md | 60 +++++++++++++++++++++-------- 1 file changed, 45 insertions(+), 15 deletions(-) diff --git a/doc/design/incoherent_ctr_models.md b/doc/design/incoherent_ctr_models.md index c9b7b20..5140bf5 100644 --- a/doc/design/incoherent_ctr_models.md +++ b/doc/design/incoherent_ctr_models.md @@ -105,14 +105,25 @@ the termination-specific replacement of exposed Film material. `F_uc` sums all of those contributions as one complex amplitude. `SXRDCrystal.F` then adds the bulk and every top-level component before the caller takes an absolute value. -For flat-height amplitudes `A_n(Q)` and exposed fractions `p_n`, the current -calculation is equivalent, up to the configured Poisson tail truncation, to +For flat-height amplitudes `A_n(Q)` and the coherent path's exposed fractions +`q_n`, the current calculation is equivalent, up to the configured Poisson +tail truncation, to ```text -A_coherent(Q) = sum_n p_n A_n(Q) +A_coherent(Q) = sum_n q_n A_n(Q) F2_coherent(Q) = abs(A_coherent(Q)) ** 2. ``` +`q_n` is written separately from the ensemble probability `p_n` used from the +next section onwards, and the two are not interchangeable. `q_n` is what +`createLayers` assigns today: `surface_occupancy` masses, upper tail folded +into the terminal bin, left unnormalized after thresholding. `p_n` is the +incoherent ensemble's `probability(n + 1)` mass over the retained interval, +renormalized to one. They agree on every interior retained bin and differ at +the boundary. That difference is exactly why the live coherent evaluation, and +not the finite state-sum reconstruction, is the authoritative `kappa = 0` +endpoint. + This is the characteristic-function form used by the classical coherent CTR roughness treatments. Harada writes the damping factor as the squared modulus of the Fourier sum over relative areas at each step height. Dale *et al.* use @@ -371,6 +382,8 @@ class IncoherentModel(LinearFitFunctions, ABC): force_recalculate=False, ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: ... def setParameters(self, values: ArrayLike) -> None: ... + def setFitParameters(self, values): # raises; use setParameters + ... def getFitErrors(self) -> np.ndarray: ... def setFitErrors(self, errors: ArrayLike | None) -> None: ... def parameter_list(self) -> list[Parameter]: ... @@ -397,6 +410,13 @@ model.addFitParameter( ) ``` +`addFitParameter` itself is the inherited `LinearFitFunctions` implementation, +which already resolves a string through `parameterLookup` and rejects a name +already present in `self.fitparnames`. Because `fitparnames` is a composite +here, that inherited check spans the wrapped crystal's names too, so the +invariant that local and coherent names stay unique after concatenation is +enforced at the point of registration rather than at `prepareFit`. + `PoissonHeightDomains` defines `parameterLookup` and `parameterLookup_inv` for its local numerical basis, just as the existing classes do. Its initial basis contains `incoherent_fraction`; calling `addFitParameter` promotes that basis @@ -732,7 +752,8 @@ fit = CTROptimizer( fit.prepareFit() assert fit.xtal is fit.model.coherent_model assert fit.n_parameters == len(fit.get_parameters()) -prepared_start = fit.startp.copy() # legacy preparation-time snapshot +prepared_start = fit.startp.copy() # legacy model-block snapshot; the full + # vector is get_parameters() assert fit.fitparnames[-(len(crystal.fitparnames) + 1)] == ( "rough_surface incoherent_fraction" ) @@ -914,8 +935,10 @@ not contain `p_n`; probability is applied only by For sorted height states, calculate the Film part cumulatively: 1. Evaluate every distinct generated Film layer amplitude once. -2. Use a prefix sum from the lowest represented height to construct the flat - Film correction for every state. +2. Use a prefix sum starting at `min(0, lowest retained state)` -- the sharp + Film boundary, or below it when the surface is etched -- to construct the + flat Film correction for every state. Starting at the lowest retained state + instead would drop the Film layers between the boundary and that state. 3. Evaluate each termination replacement once per termination-cycle state and translated height. 4. Add the appropriate termination replacement to each Film prefix. @@ -949,9 +972,11 @@ candidate support rather than from either retained mask. `PoissonHeightDomains.exact_layer_count` defaults to 10 and is a positive integer, non-fit setting. If the candidate support contains at most this many nonzero states, retain them all. For a wider support, retain a contiguous -interval containing the Poisson mode and at least `exact_layer_count` states, -then expand it using the adjacent calculated probability masses until the -profile's cumulative `tail_probability` target is met. This is a +interval containing the mode of the calculated masses -- `argmax` of +`probability(layer_numbers + 1)`, ties resolved to the lower index, never +`floor(rate)` -- and at least `exact_layer_count` states, then expand it using +the adjacent calculated probability masses until the profile's cumulative +`tail_probability` target is met. This is a probability-data cutoff, not a dependency on observed CTR values. It does not impose a ten-state cap. @@ -1012,7 +1037,11 @@ scope. 5. **Two-height anti-Bragg case:** equal populations of two height states whose amplitudes differ by a sign cancel coherently; the fully incoherent result remains the single-flat-state `F2`, and the partial model fills the - minimum linearly with `kappa`. + minimum linearly with `kappa`. This is reachable with a real profile rather + than a synthetic ensemble: `PoissonProfile(mean_change=0.5, alpha=0.0)` has + `rate == 0` and a one-half deterministic step fraction, so it populates + structural layers `-1` and `0` at exactly `0.5` each, and there + `probability(n + 1)` and `surface_occupancy(n)` coincide exactly. 6. **Nonnegative `F2`:** random valid amplitudes, probabilities, and `kappa` values never produce negative `F2` beyond roundoff. 7. **Complete-state interference:** a fixture with a nonzero bulk amplitude @@ -1196,11 +1225,12 @@ Work: Do not give the helper a single "retained mask" output. The two policies retain different sets for different reasons, and merging them is what would silently make the coherent path adopt the incoherent cutoff or the reverse. -2. Add an immutable flat-height result containing layer numbers, - normalized probabilities, raw retained mass, excluded lower/upper mass, - `iter_states()`, and an explicit diagnostic `as_array()`. Returned arrays - must be caller-owned or read-only; the production path must not require - full amplitude materialization. +2. Add an immutable flat-height result containing layer numbers, normalized + probabilities, raw retained mass, excluded lower and upper mass, their sum + as `excluded_probability`, `iter_states()`, and an explicit diagnostic + `as_array()` -- the complete attribute set sketched under "Flat-height + surface evaluation". Returned arrays must be caller-owned or read-only; the + production path must not require full amplitude materialization. 3. Add `PoissonSurface.flat_domain_corrections(h, k, l, *, exact_layer_count=10)`. For state `n`, the result is the occupancy-one component correction From 32176070b7a9b68f371b74b2dbec2183a71f2f68 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Fri, 11 Sep 2026 13:23:45 -0400 Subject: [PATCH 04/13] fix(phys): route deprecated statistics errors through set_errors evaluateStatistics hardcoded the resolution prefix and skipped every callback, so the crystal received the callback's error slice. --- orgui/datautils/xrayutils/CTRopt.py | 10 +-- orgui/datautils/xrayutils/test/test_CTRopt.py | 72 +++++++++++++++++++ 2 files changed, 77 insertions(+), 5 deletions(-) diff --git a/orgui/datautils/xrayutils/CTRopt.py b/orgui/datautils/xrayutils/CTRopt.py index fa1c781..4e7a2d7 100644 --- a/orgui/datautils/xrayutils/CTRopt.py +++ b/orgui/datautils/xrayutils/CTRopt.py @@ -1412,11 +1412,11 @@ def evaluateStatistics(self, x): pcov = util.leastsq_covariance(self.residues, x) errors = np.sqrt(np.diag(pcov) * chi2_red) - if self._fit_resolution: - self.resolution_errors = errors[:3] - self.xtal.setFitErrors(errors[3:]) - else: - self.xtal.setFitErrors(errors) + # Splitting the vector here by hand hardcoded the three-entry + # resolution prefix and skipped every registered callback, so with a + # callback present the crystal was handed the callback's error slice. + # `set_errors` is the same splitter `statistics` uses. + self.set_errors(errors) self.set_parameters(x) return chi2_result, chi2_red, pvalue, residues2.size diff --git a/orgui/datautils/xrayutils/test/test_CTRopt.py b/orgui/datautils/xrayutils/test/test_CTRopt.py index 06d7f1c..cfc6fd8 100644 --- a/orgui/datautils/xrayutils/test/test_CTRopt.py +++ b/orgui/datautils/xrayutils/test/test_CTRopt.py @@ -1636,6 +1636,78 @@ def test_set_errors_slices_by_parameter_layout(self): self.assertEqual(len(optimizer.xtal.error_calls), 1) +class TestDeprecatedStatisticsErrorRouting(unittest.TestCase): + """The legacy statistics path splits errors like ``set_errors``. + + ``evaluateStatistics`` used to call ``self.xtal.setFitErrors`` directly on + a hand-sliced vector. That hardcoded the three-entry resolution prefix and + ignored registered callbacks entirely, so the crystal was handed the + callbacks' error slice as if it were its own. + """ + + def test_callback_errors_do_not_reach_the_crystal(self): + """Each callback keeps its own slice and the crystal keeps the tail.""" + optimizer = CTRopt.CTROptimizer(FitCrystal((1.0, 2.0)), _fixture_ctrs()) + first, _ = _register_callback(optimizer, "cb_first", 2) + second, _ = _register_callback(optimizer, "cb_second", 1) + optimizer.prepareFit() + + self.assertEqual(optimizer.callbacks, [second, first]) + self.assertEqual(optimizer.get_parameters().size, 5) + + with self.assertWarnsRegex( + DeprecationWarning, "evaluateStatistics is deprecated" + ): + optimizer.evaluateStatistics(optimizer.get_parameters()) + + # The crystal tail is checked first: before the repair it received all + # five entries, so this is the assertion which names the defect. + self.assertEqual(np.size(optimizer.xtal.errors), 2) + self.assertEqual(len(optimizer.xtal.error_calls), 1) + self.assertEqual(np.size(second.errors), 1) + self.assertEqual(np.size(first.errors), 2) + + def test_resolution_prefix_is_not_hardcoded(self): + """The resolution block is split by the shared layout, not by ``[:3]``.""" + optimizer = CTRopt.CTROptimizer(FitCrystal((1.0, 2.0)), _fixture_ctrs()) + # All three widths start strictly positive. The covariance estimate + # differentiates numerically, and a width starting at zero is driven + # slightly negative by that step, which the resolution model rejects. + optimizer.fit_resolution( + CTRresolution.BoxResolution(0.1, 0.1, 0.1), + lower_bounds=[0.0, 0.0, 0.0], + higher_bounds=[1.0, 1.0, 1.0], + ) + callback, _ = _register_callback(optimizer, "cb", 2) + optimizer.prepareFit() + + self.assertEqual(optimizer.get_parameters().size, 7) + + with self.assertWarnsRegex( + DeprecationWarning, "evaluateStatistics is deprecated" + ): + optimizer.evaluateStatistics(optimizer.get_parameters()) + + # Before the repair the crystal received everything after the three + # resolution entries, the callback block included. + self.assertEqual(np.size(optimizer.xtal.errors), 2) + self.assertEqual(np.size(optimizer.resolution_errors), 3) + self.assertEqual(np.size(callback.errors), 2) + + def test_unprefixed_layout_is_unchanged(self): + """With no callbacks and no fitted resolution the crystal gets it all.""" + optimizer = CTRopt.CTROptimizer(FitCrystal((1.0, 2.0)), _fixture_ctrs()) + optimizer.prepareFit() + + with self.assertWarnsRegex( + DeprecationWarning, "evaluateStatistics is deprecated" + ): + optimizer.evaluateStatistics(optimizer.get_parameters()) + + self.assertEqual(np.size(optimizer.xtal.errors), 2) + self.assertEqual(len(optimizer.xtal.error_calls), 1) + + class TestCallbackBounds(unittest.TestCase): """Increment 1: the ``FitCallback`` bounds ordering repair.""" From ad52cbc84424fc161997cc1c7c22ef47027881df Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Fri, 11 Sep 2026 13:46:19 -0400 Subject: [PATCH 05/13] docs(phys): add commit-message cleanup as increment 7 The bodies written so far restate their own diffs. Records the subject, 0-3 line body and footer convention, to be applied by rebase while the branch is still unpushed. --- doc/design/incoherent_ctr_models.md | 37 +++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/doc/design/incoherent_ctr_models.md b/doc/design/incoherent_ctr_models.md index 5140bf5..e945bd2 100644 --- a/doc/design/incoherent_ctr_models.md +++ b/doc/design/incoherent_ctr_models.md @@ -1536,6 +1536,43 @@ Gate: - Old sessions and `.xtal` files remain untouched and load unchanged because this increment does not integrate the new dictionary contract with them. +### Increment 7 -- normalize the commit history + +Do this last, after increment 6 and before the branch is integrated or +pushed. The commits written during increments 0--6 carry multi-paragraph +bodies which restate their own diffs. They are too long. + +The convention for this repository, on top of the Conventional Commits rule +in the root `AGENTS.md`: + +- a subject line; +- optionally a body of **zero to three lines**, separated from the subject by + one blank line; +- optionally footers, separated from the body by one blank line, with tokens + using `-` in place of whitespace (`Reviewed-by:`, `Refs:`) and + `BREAKING CHANGE` as the one permitted exception; +- everything wrapped at 72 columns. + +The body says what changed and why, not how. Omit it when the change is +simple: a body which restates the subject is a sign that there should not be +one. + +Work: + +1. Reword the branch's own commits to that shape. They are unpushed, so this + is a local rebase and rewrites no published history. Check that first with + `git log origin/poisson-incoh..HEAD`; if any commit has been published, + leave it and note the exception here instead. +2. Keep every subject's existing Conventional Commits type and scope. This is + a length and layout change, not a reclassification. +3. Preserve the `BREAKING CHANGE:` footer on any commit which carries one. + +Gate: + +- No commit body exceeds three lines and no line exceeds 72 columns. +- `git log origin/poisson-incoh..HEAD` contains the same set of changes as + before the rewording, verified by comparing the tree at the branch tip. + ### Deferred increment -- physical coherence kernel Do not include this in the first implementation. A later change may add a From 5c105fd76a1bd023cbc9f9052cd64c9060db6bc1 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Fri, 11 Sep 2026 14:10:07 -0400 Subject: [PATCH 06/13] test(phys): characterize the coherent Poisson surface Records the flat-height decomposition the incoherent models are defined against, against independently stacked Films that never touch PoissonSurface. --- .../xrayutils/test/_poisson_oracle.py | 300 +++++++++++++++++ .../datautils/xrayutils/test/test_CTRcalc.py | 313 ++++++++++++++++++ 2 files changed, 613 insertions(+) create mode 100644 orgui/datautils/xrayutils/test/_poisson_oracle.py diff --git a/orgui/datautils/xrayutils/test/_poisson_oracle.py b/orgui/datautils/xrayutils/test/_poisson_oracle.py new file mode 100644 index 0000000..7498ac5 --- /dev/null +++ b/orgui/datautils/xrayutils/test/_poisson_oracle.py @@ -0,0 +1,300 @@ +"""Independent deterministic flat-height crystals for Poisson surface tests. + +Shared by ``test_CTRcalc.py`` and ``test_CTRfilm.py``. The package has no +``conftest.py`` and its CTR tests are ``unittest`` classes, so these are plain +importable builders rather than pytest fixtures. + +The oracle represents one flat surface height as a plain :class:`CTRfilm.Film` +of the corresponding thickness. It never calls ``PoissonSurface``, so a +placement or occupancy error in the Poisson assembly cannot cancel out of a +comparison against it. + +Height convention, matching ``doc/design/incoherent_ctr_models.md``: a state is +named by the structural layer ``n`` of its top filled layer. Layers ``j < 0`` +lie inside the sharp Film boundary and ``j >= 0`` above it, so state ``n`` +corresponds to a Film of ``w_base + n + 1`` layers, and the signed height +change is ``H = n + 1``. + +Units: lattice constants and heights in Angstrom, ``h``, ``k`` and ``l`` in +r.l.u. of the reference cell, amplitudes in electrons per reference lateral +cell. +""" + +import copy + +import numpy as np + +from .. import CTRcalc, CTRfilm + +__all__ = [ + "coherent_reference", + "flat_height_crystal", + "height_states", + "layered_cell", + "poisson_crystal", +] + + +def layered_cell(n_layers=2, name="layered", a=3.0, c=6.0): + """Return a unit cell whose structural layers are evenly spaced in z. + + ``UnitCell.addAtom`` defaults ``layerpos`` to ``0.0`` for every new layer, + which collapses the whole cycle onto one stacking position: every Film + layer is then placed a full unit cell apart instead of a fraction of one, + and any comparison against an independently stacked Film disagrees. The + explicit ``layerpos`` assignment below is therefore load bearing, not + decoration. + + :param int n_layers: + Number of structural layers in the cycle. + :param str name: + Unit cell name. + :param float a: + In-plane lattice constant in Angstrom. + :param float c: + Out-of-plane lattice constant in Angstrom. + :returns: + Cell with ``n_layers`` layers at fractional heights ``i / n_layers``. + :rtype: CTRcalc.UnitCell + """ + unitcell = CTRcalc.UnitCell([a, a, c], [90.0, 90.0, 90.0], name=name) + for index in range(n_layers): + z = index / n_layers + element = "C" if index % 2 == 0 else "O" + unitcell.addAtom(element, [0.0, 0.0, z], 0.1, 0.1, 1.0, layer=index) + unitcell.layerpos[float(index)] = z + return unitcell + + +def _crystal( + bulk_cell, + components, + stacking, + reference_uc=None, + weights=None, + domains=None, +): + """Assemble one crystal and apply any requested weights and domains. + + ``domains`` is applied to every top-level component, not to one of them. + A flat-height crystal has a single component where the crystal under test + has two, so applying it selectively would make the two sides differ by the + transform itself rather than by the surface model. + """ + keyargs = {"stacking": np.asarray(stacking)} + if reference_uc is not None: + keyargs["reference_uc"] = reference_uc + crystal = CTRcalc.SXRDCrystal(bulk_cell, *components, **keyargs) + if weights is not None: + crystal.weights = np.asarray(weights, dtype=np.float64) + crystal.weights_0 = np.copy(crystal.weights) + if domains is not None: + for index in range(len(components)): + crystal.setDomain(index, list(domains)) + return crystal + + +def poisson_crystal( + profile, + w_base=4.0, + n_layers=2, + termination_cells=None, + reference_uc=None, + weights=None, + domains=None, +): + """Return the crystal under test and its Poisson surface. + + :param CTRdistributions.PoissonProfile profile: + Height distribution driving the surface. + :param float w_base: + Base Film width in structural layers. + :param int n_layers: + Structural layers per unit cell. + :param dict termination_cells: + Optional explicit termination bank. Defaults to cells generated from + the Film's own cell, which makes the termination replacement cancel. + :param CTRcalc.UnitCell reference_uc: + Optional reference cell, to exercise reference-area scaling. + :param sequence weights: + Optional component weights for ``[film, surface]``. + :param sequence domains: + Optional ``(matrix, occupancy)`` list applied to every component. + :returns: + The crystal and its surface component. + :rtype: tuple[CTRcalc.SXRDCrystal, CTRfilm.PoissonSurface] + """ + film = CTRfilm.Film(layered_cell(n_layers, "film"), name="film") + film.basis[0] = w_base + source = ( + termination_cells + if termination_cells is not None + else layered_cell(n_layers, "film") + ) + surface = CTRfilm.PoissonSurface(source, profile=profile, name="surface") + crystal = _crystal( + layered_cell(n_layers, "bulk"), + (film, surface), + [1, 2], + reference_uc=reference_uc, + weights=weights, + domains=domains, + ) + return crystal, surface + + +def flat_height_crystal( + n, + w_base=4.0, + n_layers=2, + reference_uc=None, + weights=None, + domains=None, +): + """Return the independent crystal whose surface is flat at top layer ``n``. + + Built only from a plain Film of ``w_base + n + 1`` layers, with no + ``PoissonSurface`` anywhere. + + :param int n: + Structural layer of the top filled layer. + :param float w_base: + Base Film width in structural layers. + :param int n_layers: + Structural layers per unit cell. + :param CTRcalc.UnitCell reference_uc: + Optional reference cell, matching ``poisson_crystal``. + :param sequence weights: + Optional single-element component weight for the Film. + :param sequence domains: + Optional ``(matrix, occupancy)`` list applied to the Film. + :returns: + Crystal with a sharp surface at the requested height. + :rtype: CTRcalc.SXRDCrystal + :raises ValueError: + If the requested height etches away more than the base Film. + """ + width = w_base + n + 1 + if width < 0: + raise ValueError( + f"flat height {n} etches below the base film of {w_base} layers" + ) + film = CTRfilm.Film(layered_cell(n_layers, "film"), name="film") + film.basis[0] = width + return _crystal( + layered_cell(n_layers, "bulk"), + (film,), + [1], + reference_uc=reference_uc, + weights=weights, + domains=domains, + ) + + +def height_states(profile, minimum_probability=0.0): + """Return the represented structural layers and their exposed fractions. + + These are the coherent path's ``q_n``: ``surface_occupancy`` masses, with + the upper tail folded into the terminal bin. They are not the normalized + ensemble probabilities of the incoherent model. + + :param CTRdistributions.PoissonProfile profile: + Height distribution. + :param float minimum_probability: + Discard states at or below this exposed fraction. + :returns: + Structural layer numbers and their exposed fractions. + :rtype: tuple[numpy.ndarray, numpy.ndarray] + """ + lower, upper = profile.support() + layers = np.arange(lower, upper + 1) + exposed = profile.surface_occupancy(layers) + keep = exposed > minimum_probability + return layers[keep], exposed[keep] + + +def coherent_reference(profile, h, k, l, **keyargs): # noqa: E741 + """Return ``sum_n q_n A_n`` over independent flat-height crystals. + + :param CTRdistributions.PoissonProfile profile: + Height distribution. + :param numpy.ndarray h: + Reference-frame reciprocal coordinate in r.l.u. + :param numpy.ndarray k: + Reference-frame reciprocal coordinate in r.l.u. + :param numpy.ndarray l: + Reference-frame reciprocal coordinate in r.l.u. + :param keyargs: + Forwarded to :func:`flat_height_crystal`. + :returns: + Coherent amplitude average in electrons per reference lateral cell. + :rtype: numpy.ndarray + """ + layers, exposed = height_states(profile) + total = np.zeros_like(np.asarray(l, dtype=np.float64), dtype=np.complex128) + for layer, probability in zip(layers, exposed): + crystal = flat_height_crystal(int(layer), **keyargs) + total = total + probability * crystal.F(h, k, l) + return total + + +def flat_height_amplitudes(profile, h, k, l, **keyargs): # noqa: E741 + """Return the layers, exposed fractions, and per-state amplitudes. + + The squared-amplitude average ``sum_n q_n abs(A_n) ** 2`` built from these + is the fully incoherent limit which later increments must reproduce. + + :returns: + Structural layers, exposed fractions, and one amplitude array per + state. + :rtype: tuple[numpy.ndarray, numpy.ndarray, list[numpy.ndarray]] + """ + layers, exposed = height_states(profile) + amplitudes = [ + flat_height_crystal(int(layer), **keyargs).F(h, k, l) + for layer in layers + ] + return layers, exposed, amplitudes + + +def minimum_base_width(profile, margin=2.0): + """Return a base Film width which survives the deepest etched state. + + :param CTRdistributions.PoissonProfile profile: + Height distribution. + :param float margin: + Extra layers kept below the deepest state. + :returns: + Base width in structural layers. + :rtype: float + """ + layers, _ = height_states(profile) + return float(max(margin, -int(layers[0]) + margin)) + + +def termination_bank(n_layers=2, scale=1.0, displacement=0.0): + """Return a distinct termination bank for every layer of the cycle. + + :param int n_layers: + Structural layers per unit cell. + :param float scale: + Occupancy multiplier applied to every termination atom. + :param float displacement: + Fractional z displacement applied to every termination atom. + :returns: + Mapping from structural layer to termination cell. + :rtype: dict + """ + source = layered_cell(n_layers, "termination") + if scale != 1.0 or displacement != 0.0: + source.basis[:, 6] = source.basis[:, 6] * scale + source.basis[:, 3] = source.basis[:, 3] + displacement + film_cell = layered_cell(n_layers, "film") + return CTRfilm.generate_surface_termination_cells( + source, np.asarray(list(film_cell.split_in_layers())) + ) + + +def deep_copy_cells(cells): + """Return an independent copy of a termination bank.""" + return {layer: copy.deepcopy(cell) for layer, cell in cells.items()} diff --git a/orgui/datautils/xrayutils/test/test_CTRcalc.py b/orgui/datautils/xrayutils/test/test_CTRcalc.py index 44c8fa9..c22e1a0 100644 --- a/orgui/datautils/xrayutils/test/test_CTRcalc.py +++ b/orgui/datautils/xrayutils/test/test_CTRcalc.py @@ -44,6 +44,7 @@ from .. import CTRcalc, CTRfilm, CTRplotutil, CTRsymmetry, CTRuc from ..CTRdistributions import PoissonProfile, SkellamProfile, SurfaceProfile from ..CTRutil import generate_surface_termination_cells +from . import _poisson_oracle HAS_ASE = importlib.util.find_spec("ase") is not None @@ -1039,6 +1040,318 @@ def test_repeated_stacking_refreshes_underlying_film_views(self): 0.5, ) + +class TestPoissonFlatHeightCharacterization(unittest.TestCase): + """Characterize the coherent Poisson surface against flat-height crystals. + + The reference in every test is an independently stacked ``Film`` of the + matching thickness, built by ``_poisson_oracle`` without touching + ``PoissonSurface``. These record the behavior the incoherent models of + ``doc/design/incoherent_ctr_models.md`` are defined against; they are + characterization tests, so a deliberate change to the coherent model is + expected to update them. + """ + + H = np.zeros(5) + K = np.zeros(5) + L = np.array([0.35, 0.8, 1.3, 1.85, 2.4]) + + def relative_error(self, value, reference): + """Return the max deviation relative to the reference amplitude.""" + return float( + np.max(np.abs(value - reference)) / np.max(np.abs(reference)) + ) + + def test_layered_cell_fixture_sets_every_layer_position(self): + """``addAtom`` defaults ``layerpos`` to zero for each new layer. + + A cell left that way collapses its whole cycle onto one stacking + position, which places consecutive Film layers a full unit cell apart. + Every comparison in this class then disagrees with the Poisson surface + for reasons that have nothing to do with the surface, so the fixture + pins the positions explicitly. + """ + collapsed = CTRcalc.UnitCell([3.0, 3.0, 6.0], [90.0, 90.0, 90.0]) + for index, z in enumerate((0.0, 0.5)): + collapsed.addAtom("C", [0.0, 0.0, z], 0.1, 0.1, 1.0, layer=index) + self.assertEqual(collapsed.layerpos, {0.0: 0.0, 1.0: 0.0}) + + self.assertEqual( + _poisson_oracle.layered_cell(2).layerpos, + {0.0: 0.0, 1.0: 0.5}, + ) + + def test_deterministic_states_match_independent_films(self): + """One flat height equals a plain Film of the matching thickness. + + Covers positive growth and negative etching, layer cycles which do and + do not divide the base width, and cycles of one, two, and three + structural layers. + """ + for n_layers in (1, 2, 3): + for w_base in (3.0, 4.0, 5.0): + for height_change in (2, 1, 0, -1, -2): + with self.subTest( + n_layers=n_layers, + w_base=w_base, + height_change=height_change, + ): + profile = PoissonProfile( + float(height_change), + alpha=0.0, + offset=0.0, + tail_probability=1e-14, + ) + crystal, _ = _poisson_oracle.poisson_crystal( + profile, w_base=w_base, n_layers=n_layers + ) + reference = _poisson_oracle.flat_height_crystal( + height_change - 1, + w_base=w_base, + n_layers=n_layers, + ) + self.assertLess( + self.relative_error( + crystal.F(self.H, self.K, self.L), + reference.F(self.H, self.K, self.L), + ), + 1e-12, + ) + + def test_coherent_amplitude_is_the_flat_height_average(self): + """``F`` equals ``sum_n q_n A_n`` over independent flat crystals. + + This is the characteristic-function form the design record calls the + coherent limit. It is recorded for growth and etching, for a + fractional deterministic step, for a nonzero offset, and for mixtures + of a step with a Poisson tail. + """ + cases = { + "integer growth": PoissonProfile(2.0, alpha=0.0), + "fractional step": PoissonProfile(1.5, alpha=0.0), + "nonzero offset": PoissonProfile(1.0, alpha=0.0, offset=0.25), + "negative offset": PoissonProfile(0.0, alpha=0.0, offset=-2.25), + "poisson growth": PoissonProfile(1.5, alpha=0.6), + "poisson etching": PoissonProfile(-1.5, alpha=0.6), + "step and poisson": PoissonProfile(2.3, alpha=0.4, offset=0.4), + "deep etching": PoissonProfile(-4.0, alpha=1.0), + } + for label, profile in cases.items(): + with self.subTest(case=label): + profile.tail_probability = 1e-14 + w_base = _poisson_oracle.minimum_base_width(profile) + crystal, _ = _poisson_oracle.poisson_crystal( + profile, w_base=w_base + ) + reference = _poisson_oracle.coherent_reference( + profile, self.H, self.K, self.L, w_base=w_base + ) + self.assertLess( + self.relative_error( + crystal.F(self.H, self.K, self.L), reference + ), + 1e-12, + ) + + def test_incoherent_average_departs_from_the_coherent_limit(self): + """Two equally populated heights cancel coherently but not squared. + + ``PoissonProfile(0.5, alpha=0.0)`` has zero Poisson rate and a + one-half deterministic step fraction, so it populates structural + layers ``-1`` and ``0`` at exactly one half each. The coherent square + of the averaged amplitude and the averaged square of the two state + amplitudes are the two endpoints the incoherent model interpolates, + and they must differ here. + """ + profile = PoissonProfile(0.5, alpha=0.0, tail_probability=1e-14) + layers, exposed = _poisson_oracle.height_states(profile) + np.testing.assert_array_equal(layers, [-1, 0]) + np.testing.assert_allclose(exposed, [0.5, 0.5]) + + # The contrast between the two endpoints varies strongly along the + # rod, so it is scanned rather than sampled at a few points: the + # bulk dominates near the Bragg conditions and suppresses it there. + scan = np.linspace(0.05, 3.0, 296) + zeros = np.zeros_like(scan) + _, exposed, amplitudes = _poisson_oracle.flat_height_amplitudes( + profile, zeros, zeros, scan, w_base=4.0 + ) + coherent = sum( + probability * amplitude + for probability, amplitude in zip(exposed, amplitudes) + ) + incoherent = sum( + probability * np.abs(amplitude) ** 2 + for probability, amplitude in zip(exposed, amplitudes) + ) + + crystal, _ = _poisson_oracle.poisson_crystal(profile, w_base=4.0) + np.testing.assert_allclose( + crystal.F(zeros, zeros, scan), coherent, rtol=1e-12 + ) + + # The incoherent endpoint is never below the coherent one: it drops + # the cross term rather than adding anything, so their difference is + # the state variance. + variance = incoherent - np.abs(coherent) ** 2 + self.assertTrue(np.all(variance >= -1e-9)) + self.assertGreater(np.max(variance / incoherent), 0.15) + + # A distribution with one populated height has no variance to expose, + # which is the invariance the partial model must preserve for every + # mixing fraction. + single = PoissonProfile(1.0, alpha=0.0, tail_probability=1e-14) + _, single_exposed, single_amplitudes = ( + _poisson_oracle.flat_height_amplitudes( + single, zeros, zeros, scan, w_base=4.0 + ) + ) + self.assertEqual(len(single_amplitudes), 1) + np.testing.assert_allclose( + single_exposed[0] * np.abs(single_amplitudes[0]) ** 2, + np.abs(single_exposed[0] * single_amplitudes[0]) ** 2, + rtol=1e-12, + ) + + def test_state_squares_keep_the_bulk_inside_each_state(self): + """Averaging only the surface correction is a different quantity. + + A fixture with a nonzero bulk amplitude separates the correct + ``sum_n q_n abs(A_n) ** 2`` from the incorrect + ``abs(A_common) ** 2 + sum_n q_n abs(dA_n) ** 2``, which drops the + bulk-surface and Film-surface interference inside each domain. + """ + profile = PoissonProfile(0.5, alpha=0.0, tail_probability=1e-14) + w_base = 4.0 + _, exposed, amplitudes = _poisson_oracle.flat_height_amplitudes( + profile, self.H, self.K, self.L, w_base=w_base + ) + common = _poisson_oracle.flat_height_crystal( + -1, w_base=w_base + ).F(self.H, self.K, self.L) + + correct = sum( + probability * np.abs(amplitude) ** 2 + for probability, amplitude in zip(exposed, amplitudes) + ) + incorrect = np.abs(common) ** 2 + sum( + probability * np.abs(amplitude - common) ** 2 + for probability, amplitude in zip(exposed, amplitudes) + ) + + self.assertGreater(np.max(np.abs(common)), 1.0) + self.assertGreater( + np.max(np.abs(correct - incorrect)) / np.max(correct), 0.1 + ) + + def test_crystal_level_scaling_preserves_the_flat_height_average(self): + """Area scaling, weights, and domain transforms stay outside the sum. + + The design record requires each state amplitude to carry the crystal's + reference-area scaling, component weight, and outer coherent-domain + transforms. Each is applied to both sides here, so the flat-height + average must survive all three and their combination. + """ + profile = PoissonProfile(1.5, alpha=0.5, tail_probability=1e-14) + w_base = _poisson_oracle.minimum_base_width(profile) + wide_reference = _poisson_oracle.layered_cell(2, "reference", a=6.0) + transforms = [ + (np.identity(3), 0.65), + (np.diag([1.0, 1.0, 0.5]), 0.35), + ] + cases = { + "reference area": {"reference_uc": wide_reference}, + "component weight": {"weights": 0.4}, + "domain transform": {"domains": transforms}, + "combined": { + "reference_uc": wide_reference, + "weights": 0.4, + "domains": transforms, + }, + } + plain, _ = _poisson_oracle.poisson_crystal(profile, w_base=w_base) + unscaled = plain.F(self.H, self.K, self.L) + for label, keyargs in cases.items(): + with self.subTest(case=label): + weight = keyargs.pop("weights", None) + crystal, _ = _poisson_oracle.poisson_crystal( + profile, + w_base=w_base, + # One weight per component: the Film and the surface + # correction together make up one flat-height Film. + weights=None if weight is None else [weight, weight], + **keyargs, + ) + reference = _poisson_oracle.coherent_reference( + profile, + self.H, + self.K, + self.L, + w_base=w_base, + weights=None if weight is None else [weight], + **keyargs, + ) + if weight is not None: + keyargs["weights"] = weight + scaled = crystal.F(self.H, self.K, self.L) + self.assertLess( + self.relative_error(scaled, reference), 1e-12 + ) + # Each factor must actually move the amplitude, or the + # agreement above would hold for a crystal ignoring it. + self.assertGreater( + self.relative_error(scaled, unscaled), 0.1 + ) + + def test_distinct_terminations_stay_linear_in_the_mixture(self): + """A termination bank does not break linearity in the height states. + + The flat-height Film oracle cannot represent a termination which + differs from the Film's own cell, so this compares the mixture against + single-state Poisson surfaces instead. It therefore characterizes + linearity rather than absolute placement, which is what the state-sum + decomposition of the incoherent model relies on. + """ + bank = _poisson_oracle.termination_bank( + 2, scale=0.75, displacement=0.05 + ) + mixed = PoissonProfile(0.5, alpha=0.0, tail_probability=1e-14) + layers, exposed = _poisson_oracle.height_states(mixed) + np.testing.assert_array_equal(layers, [-1, 0]) + + blended = np.zeros_like(self.L, dtype=np.complex128) + for layer, probability in zip(layers, exposed): + single = PoissonProfile( + float(layer) + 1.0, alpha=0.0, tail_probability=1e-14 + ) + crystal, _ = _poisson_oracle.poisson_crystal( + single, + w_base=4.0, + termination_cells=_poisson_oracle.deep_copy_cells(bank), + ) + blended = blended + probability * crystal.F( + self.H, self.K, self.L + ) + + crystal, surface = _poisson_oracle.poisson_crystal( + mixed, + w_base=4.0, + termination_cells=_poisson_oracle.deep_copy_cells(bank), + ) + mixture = crystal.F(self.H, self.K, self.L) + self.assertLess(self.relative_error(mixture, blended), 1e-12) + + # The bank must actually change the result, or the comparison above + # would pass against a surface which ignored it entirely. + plain, _ = _poisson_oracle.poisson_crystal(mixed, w_base=4.0) + self.assertGreater( + self.relative_error( + mixture, plain.F(self.H, self.K, self.L) + ), + 1e-3, + ) + + class TestLayerStacking(unittest.TestCase): @staticmethod def make_layered_unitcell(name="layered"): From 27d297987ceeb6d6770230624eb6ac70dac28b9b Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Fri, 11 Sep 2026 17:45:14 -0400 Subject: [PATCH 07/13] feat(phys): add a kinematical F2 boundary and decomposition F is now a thin wrapper over evaluate_kinematic, so incoherent models reuse the shared bulk and Film amplitudes instead of re-evaluating them. Verified bitwise identical to the previous F over 58 fixtures. --- orgui/datautils/xrayutils/CTRcalc.py | 192 +++++++++++++++++- .../datautils/xrayutils/test/test_CTRcalc.py | 121 +++++++++++ 2 files changed, 304 insertions(+), 9 deletions(-) diff --git a/orgui/datautils/xrayutils/CTRcalc.py b/orgui/datautils/xrayutils/CTRcalc.py index 38ab076..903819c 100644 --- a/orgui/datautils/xrayutils/CTRcalc.py +++ b/orgui/datautils/xrayutils/CTRcalc.py @@ -32,6 +32,7 @@ from .. import util import warnings import os +from dataclasses import dataclass # random.seed(45) import errno @@ -58,6 +59,81 @@ ) +@dataclass(frozen=True) +class KinematicComponentAmplitude: + """One top-level component's scaled contribution to ``SXRDCrystal.F``. + + :param int index: + Position in ``SXRDCrystal.uc_surface_list``. + :param str name: + Stable component name. Component selection uses this rather than + object identity, which does not survive ``copy.deepcopy``. + :param amplitude: + Complex amplitude in electrons per reference lateral cell, already + carrying the reference-area scaling, the component weight, and every + outer coherent-domain transform and occupancy. + """ + + index: int + name: str + amplitude: object + + +@dataclass(frozen=True) +class KinematicAmplitudeResult: + """Immutable decomposition of one kinematical ``SXRDCrystal`` evaluation. + + ``total`` is accumulated contribution by contribution in crystal order, so + it reproduces ``SXRDCrystal.F`` bitwise. It therefore agrees with + ``bulk + sum(part.amplitude for part in components)`` to floating-point + rounding rather than exactly. + + The arrays are caller-owned: mutating one does not affect the crystal. + When the crystal has no surface components, ``total`` and ``bulk`` are the + same array. + + :param bulk: + Semi-infinite bulk amplitude including reference-area scaling and + attenuation, in electrons per reference lateral cell. + :param tuple components: + One :class:`KinematicComponentAmplitude` per top-level component, in + ``SXRDCrystal.uc_surface_list`` order. + :param total: + Complete crystal amplitude, identical to ``SXRDCrystal.F``. + """ + + bulk: object + components: tuple + total: object + + def component(self, name): + """Return the amplitude record of the component with this name. + + :param str name: + Stable component name. + :returns: + The matching component record. + :rtype: KinematicComponentAmplitude + :raises KeyError: + If no component carries that name. + :raises ValueError: + If more than one does, which would make selection by name + ambiguous. + """ + matches = [part for part in self.components if part.name == name] + if not matches: + available = ", ".join(part.name for part in self.components) + raise KeyError( + f"No component named {name!r}. Available: {available}" + ) + if len(matches) > 1: + raise ValueError( + f"{len(matches)} components are named {name!r}; component " + "names must be unique to select one" + ) + return matches[0] + + class SXRDCrystal: """Compose bulk and surface amplitudes on one reference lateral cell. @@ -325,20 +401,118 @@ def F(self, harray, karray, Larray): Complex crystal amplitude in electrons per reference lateral cell. :rtype: numpy.ndarray """ + return self.evaluate_kinematic(harray, karray, Larray).total + + def F2(self, harray, karray, Larray): + """Return the squared crystal structure factor. + + This is ``abs(F) ** 2`` in the squared units of :meth:`F`, nominally + electrons squared per reference lateral cell. It is not detector + counts: incident flux, footprint, polarization and Lorentz factors, + detector response, acquisition time, background, and any fitted + experimental scale are all outside it. It is also not a reflectivity, + which is a dimensionless intensity ratio formed from an optical + reflection amplitude. + + :param numpy.ndarray harray: + Reference-frame reciprocal coordinate in r.l.u. + :param numpy.ndarray karray: + Reference-frame reciprocal coordinate in r.l.u. + :param numpy.ndarray Larray: + Reference-frame reciprocal coordinate in r.l.u. + :returns: + Real, nonnegative squared structure factor. + :rtype: numpy.ndarray + """ + return np.abs(self.F(harray, karray, Larray)) ** 2 + + @staticmethod + def _component_F_uc(uc, harray, karray, Larray): + """Evaluate one component's own amplitude at the given coordinates.""" + return uc.F_uc(harray, karray, Larray) + + def _iter_scaled_contributions(self, index, hkl, evaluator): + """Yield one scaled contribution per outer coherent domain. + + Each domain matrix transforms the coordinates *before* ``evaluator`` + is called, so a caller evaluating something other than the component's + own ``F_uc`` -- one flat height state, say -- still sees the + coordinates that domain requires. The yielded values already carry the + reference-area scaling, the component weight, and the domain + occupancy. + + :param int index: + Position in :attr:`uc_surface_list`. + :param numpy.ndarray hkl: + Stacked ``(3, N)`` reference-frame coordinates in r.l.u. + :param callable evaluator: + Called as ``evaluator(uc, h, k, l)`` and returning an amplitude in + electrons per the component's own lateral cell. + """ + uc = self.uc_surface_list[index] + area_scale = self.reference_area / uc.uc_area + weight = self.weights[index] + for matrix, occup in self.domains[index]: + hkl_n = np.dot(matrix, hkl) + yield area_scale * occup * weight * evaluator( + uc, hkl_n[0], hkl_n[1], hkl_n[2] + ) + + def evaluate_kinematic(self, harray, karray, Larray): + """Return the decomposed kinematical evaluation behind :meth:`F`. + + Coordinates and units are exactly those of :meth:`F`. The bulk + amplitude carries the reference-area scaling, and every component + amplitude carries its own area scaling, weight, and outer + coherent-domain transforms and occupancies. + + :param numpy.ndarray harray: + Reference-frame reciprocal coordinate in r.l.u. + :param numpy.ndarray karray: + Reference-frame reciprocal coordinate in r.l.u. + :param numpy.ndarray Larray: + Reference-frame reciprocal coordinate in r.l.u. + :returns: + Bulk amplitude, per-component amplitudes in crystal order, and + their total. + :rtype: KinematicAmplitudeResult + """ bulk_scale = self.reference_area / self.uc_bulk.uc_area - F = bulk_scale * self.uc_bulk.F_bulk(harray, karray, Larray, self.atten) + bulk = bulk_scale * self.uc_bulk.F_bulk( + harray, karray, Larray, self.atten + ) hkl = np.vstack((harray, karray, Larray)) if self.enable_uc_stacking: self.apply_stacking() - for uc, weight, domains in zip( - self.uc_surface_list, self.weights, self.domains - ): - area_scale = self.reference_area / uc.uc_area - for matrix, occup in domains: - hkl_n = np.dot(matrix, hkl) - F += area_scale * occup * weight * uc.F_uc(hkl_n[0], hkl_n[1], hkl_n[2]) - return F + # `total` is accumulated contribution by contribution in the original + # component-then-domain order, so it reproduces the pre-refactor `F` + # bitwise rather than only to rounding. + total = bulk + components = [] + for index, uc in enumerate(self.uc_surface_list): + amplitude = None + for contribution in self._iter_scaled_contributions( + index, hkl, self._component_F_uc + ): + total = total + contribution + amplitude = ( + contribution + if amplitude is None + else amplitude + contribution + ) + if amplitude is None: + amplitude = np.zeros_like( + np.asarray(bulk), dtype=np.complex128 + ) + components.append( + KinematicComponentAmplitude( + index, + getattr(uc, "name", f"component_{index}"), + amplitude, + ) + ) + return KinematicAmplitudeResult(bulk, tuple(components), total) def setDomain(self, uc_no, domains): """ diff --git a/orgui/datautils/xrayutils/test/test_CTRcalc.py b/orgui/datautils/xrayutils/test/test_CTRcalc.py index c22e1a0..97ababf 100644 --- a/orgui/datautils/xrayutils/test/test_CTRcalc.py +++ b/orgui/datautils/xrayutils/test/test_CTRcalc.py @@ -1352,6 +1352,127 @@ def test_distinct_terminations_stay_linear_in_the_mixture(self): ) +class TestKinematicDecomposition(unittest.TestCase): + """``evaluate_kinematic`` exposes the pieces ``F`` already summed.""" + + H = np.zeros(5) + K = np.zeros(5) + L = np.array([0.35, 0.8, 1.3, 1.85, 2.4]) + + def crystal(self, **keyargs): + """Return a two-component crystal with a rough surface.""" + profile = PoissonProfile(1.5, alpha=0.5, tail_probability=1e-14) + crystal, _ = _poisson_oracle.poisson_crystal( + profile, w_base=4.0, **keyargs + ) + return crystal + + def test_total_reproduces_F_exactly(self): + """``total`` is the same accumulation ``F`` performed before.""" + crystal = self.crystal() + result = crystal.evaluate_kinematic(self.H, self.K, self.L) + np.testing.assert_array_equal( + result.total, crystal.F(self.H, self.K, self.L) + ) + + def test_parts_sum_to_the_total(self): + """No contribution is missing from the decomposition.""" + crystal = self.crystal() + result = crystal.evaluate_kinematic(self.H, self.K, self.L) + rebuilt = result.bulk + sum( + part.amplitude for part in result.components + ) + # Equal to rounding, not bitwise: `total` keeps the historical + # accumulation order, which interleaves the components differently. + np.testing.assert_allclose(result.total, rebuilt, rtol=1e-14) + + def test_bulk_matches_a_crystal_without_components(self): + """The bulk term carries its reference-area scaling and attenuation.""" + crystal = self.crystal() + result = crystal.evaluate_kinematic(self.H, self.K, self.L) + bulk_only = CTRcalc.SXRDCrystal(_poisson_oracle.layered_cell(2, "bulk")) + np.testing.assert_allclose( + result.bulk, + bulk_only.F(self.H, self.K, self.L), + rtol=1e-14, + ) + + def test_components_keep_crystal_order_and_names(self): + """Components are reported in ``uc_surface_list`` order.""" + crystal = self.crystal() + result = crystal.evaluate_kinematic(self.H, self.K, self.L) + self.assertEqual( + [(part.index, part.name) for part in result.components], + [(0, "film"), (1, "surface")], + ) + self.assertIs(result.component("surface"), result.components[1]) + + def test_component_lookup_rejects_unknown_and_ambiguous_names(self): + """Selection is by stable name, so the name must identify one part.""" + crystal = self.crystal() + result = crystal.evaluate_kinematic(self.H, self.K, self.L) + with self.assertRaisesRegex(KeyError, "No component named"): + result.component("missing") + + crystal.uc_surface_list[1].name = "film" + clashing = crystal.evaluate_kinematic(self.H, self.K, self.L) + with self.assertRaisesRegex(ValueError, "must be unique"): + clashing.component("film") + + def test_component_amplitudes_carry_weight_and_domains(self): + """Weights and outer domain transforms are inside each component.""" + plain = self.crystal().evaluate_kinematic(self.H, self.K, self.L) + weighted = self.crystal( + weights=[1.0, 0.25] + ).evaluate_kinematic(self.H, self.K, self.L) + + np.testing.assert_allclose( + weighted.component("surface").amplitude, + 0.25 * plain.component("surface").amplitude, + rtol=1e-14, + ) + np.testing.assert_allclose( + weighted.component("film").amplitude, + plain.component("film").amplitude, + rtol=1e-14, + ) + + halved = self.crystal( + domains=[(np.identity(3), 0.5)] + ).evaluate_kinematic(self.H, self.K, self.L) + np.testing.assert_allclose( + halved.component("surface").amplitude, + 0.5 * plain.component("surface").amplitude, + rtol=1e-14, + ) + + def test_F2_is_the_squared_modulus_of_F(self): + """``F2`` adds a squared boundary without changing the quantity.""" + crystal = self.crystal() + amplitude = crystal.F(self.H, self.K, self.L) + squared = crystal.F2(self.H, self.K, self.L) + np.testing.assert_allclose( + squared, np.abs(amplitude) ** 2, rtol=1e-14 + ) + self.assertTrue(np.all(np.isreal(squared))) + self.assertTrue(np.all(squared >= 0.0)) + + # Scalar coordinates are rejected inside `F_bulk`, which predates this + # boundary; `F2` must inherit that rather than diverge from `F`. + with self.assertRaises(AttributeError): + crystal.F(0.0, 0.0, 1.3) + with self.assertRaises(AttributeError): + crystal.F2(0.0, 0.0, 1.3) + + one_point = np.array([1.3]) + np.testing.assert_allclose( + crystal.F2(one_point * 0.0, one_point * 0.0, one_point), + np.abs(crystal.F(one_point * 0.0, one_point * 0.0, one_point)) + ** 2, + rtol=1e-14, + ) + + class TestLayerStacking(unittest.TestCase): @staticmethod def make_layered_unitcell(name="layered"): From 25b3636197f406f307277618771e4c1b5d479ca9 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Fri, 11 Sep 2026 18:39:27 -0400 Subject: [PATCH 08/13] feat(phys): decompose the Poisson surface into flat heights flat_domain_corrections returns the occupancy-one correction of each retained height, streaming one amplitude at a time. Coherent assembly and the height-state ensemble now select from shared candidates. --- orgui/datautils/xrayutils/CTRfilm.py | 625 ++++++++++++++++-- .../datautils/xrayutils/test/test_CTRfilm.py | 354 ++++++++++ 2 files changed, 907 insertions(+), 72 deletions(-) create mode 100644 orgui/datautils/xrayutils/test/test_CTRfilm.py diff --git a/orgui/datautils/xrayutils/CTRfilm.py b/orgui/datautils/xrayutils/CTRfilm.py index 13ee561..66e489a 100644 --- a/orgui/datautils/xrayutils/CTRfilm.py +++ b/orgui/datautils/xrayutils/CTRfilm.py @@ -32,6 +32,7 @@ from .. import util import re from collections.abc import Mapping +from dataclasses import dataclass, field # random.seed(45) @@ -83,6 +84,180 @@ def _unwrapped_layer_positions(layer_positions, order): return positions +@dataclass(frozen=True) +class _PoissonCandidates: + """Unmasked profile arrays for one Poisson surface evaluation. + + No retention policy is applied here. The coherent layer assembly and the + incoherent height-state ensemble select from these arrays separately, + because they retain different sets for different reasons: the coherent + path keeps every layer carrying either exposed surface or a Film + correction and never renormalizes, while the ensemble keeps exposed + heights only and renormalizes over them. + + :param PoissonProfile profile: + Profile rebuilt from the live basis. + :param numpy.ndarray layer_numbers: + Consecutive structural layer offsets covering the profile support. + :param numpy.ndarray material_occupancy: + Cumulative material occupancy, ``P(H > n)`` for height change ``H``. + :param numpy.ndarray surface_occupancy: + Exposed fraction per layer as the coherent assembly uses it. Its + terminal bin carries the folded upper tail. + :param numpy.ndarray exposed_probability: + Height-state mass ``probability(n + 1)``. Layer ``n`` is the top + filled layer exactly when ``H == n + 1``, which is why the argument is + shifted. This agrees with ``surface_occupancy`` on every interior bin + and deliberately differs at the terminal one. + :param numpy.ndarray film_correction_occupancy: + Material occupancy relative to the sharp Film boundary. + :param float tail_probability: + Cumulative tail target carried by the profile. + """ + + profile: object + layer_numbers: np.ndarray + material_occupancy: np.ndarray + surface_occupancy: np.ndarray + exposed_probability: np.ndarray + film_correction_occupancy: np.ndarray + tail_probability: float + + +@dataclass(frozen=True) +class FlatHeightState: + """One flat surface height and its occupancy-one correction. + + :param int layer_number: + Structural layer of the top filled layer. + :param float probability: + Normalized mass of this state within the retained interval. + :param amplitude: + Complex correction relative to the sharp Film boundary, in electrons + for one lateral surface unit cell. The probability is *not* folded + into it. + """ + + layer_number: int + probability: float + amplitude: object + + +@dataclass(frozen=True) +class FlatHeightCorrections: + """Retained flat-height states of a Poisson surface. + + ``probabilities`` are renormalized over the retained interval exactly + once, after selection, so they sum to one. The masses excluded below and + above the interval are reported separately and are taken from the + profile's own cumulative distribution, so the three add to one regardless + of where the candidate support was truncated. + + The stored arrays are read-only. Amplitudes are produced lazily by + :meth:`iter_states` and are only valid while the surface's parameters are + unchanged. + + :param numpy.ndarray layer_numbers: + Retained structural layers in ascending order. + :param numpy.ndarray probabilities: + Normalized state masses aligned with ``layer_numbers``. + :param float raw_retained_probability: + Retained mass before normalization. + :param float excluded_lower_probability: + Mass of heights below the retained interval. + :param float excluded_upper_probability: + Mass of heights above the retained interval. + :param callable amplitude_source: + Internal. Called with no arguments to obtain a generator of + ``(layer_number, amplitude)`` pairs. + """ + + layer_numbers: np.ndarray + probabilities: np.ndarray + raw_retained_probability: float + excluded_lower_probability: float + excluded_upper_probability: float + amplitude_source: object = field(repr=False) + + def __post_init__(self): + for name in ("layer_numbers", "probabilities"): + array = getattr(self, name) + array.flags.writeable = False + + @property + def excluded_probability(self): + """Return the total mass outside the retained interval.""" + return ( + self.excluded_lower_probability + self.excluded_upper_probability + ) + + def iter_states(self): + """Yield one :class:`FlatHeightState` at a time. + + States are produced outward from the sharp Film boundary: those at or + above it in ascending order, then those below it in descending order. + That order lets the cumulative Film correction be built with one + evaluation per structural layer and one amplitude array in flight, so + it is the evaluation order rather than the ascending order of + :attr:`layer_numbers`. Each state carries its own layer number and + probability, so a consumer accumulating sums does not depend on it. + + :returns: + Generator of flat-height states. + """ + masses = dict(zip(self.layer_numbers.tolist(), self.probabilities)) + for layer_number, amplitude in self.amplitude_source(): + yield FlatHeightState( + int(layer_number), float(masses[layer_number]), amplitude + ) + + def as_array(self): + """Return every state amplitude, in ``layer_numbers`` order. + + Diagnostic only: this materializes an ``(n_states, n_points)`` complex + array which the production path deliberately avoids. + + :returns: + Stacked state corrections aligned with :attr:`layer_numbers`. + :rtype: numpy.ndarray + """ + collected = { + int(state.layer_number): state.amplitude + for state in self.iter_states() + } + return np.array( + [collected[int(layer)] for layer in self.layer_numbers] + ) + + +def _retained_mass(candidates, low, high): + """Return the exact mass of the retained height-state interval. + + Taken from the profile's cumulative distribution rather than by summing + the candidate masses, so truncation of the candidate support does not + leak into the reported retained and excluded masses. + """ + layers = candidates.layer_numbers + occupancy = candidates.profile.occupancy( + [layers[low], layers[high] + 1] + ) + return float(occupancy[0] - occupancy[1]) + + +@dataclass(frozen=True) +class _LayerGeometry: + """Placement of the three slabs belonging to one structural layer.""" + + layer_number: int + layer_id: float + surface_uc: object + surface_matrix: np.ndarray + film_uc: object + film_matrix: np.ndarray + reference_uc: object + reference_matrix: np.ndarray + + def _translate_domains(layers, height): for layer in layers: offset = height / layer.a[2] @@ -2233,6 +2408,366 @@ def mean_height_absolute(self): layer_height = film_uc.a[2] / len(self._layer_ids) return self.below_H + layer_height * (self.basis[0] + self.basis[2]) + def _profile_candidates(self): + """Return the unmasked profile arrays for the live basis. + + Applies no retention policy. The coherent assembly and the incoherent + height-state ensemble select from the result separately. + + :returns: + Candidate layer numbers, occupancies, and masses. + :rtype: _PoissonCandidates + """ + tail_probability = ( + self.profile.tail_probability + if self.profile is not None + else DEFAULT_TAIL_PROBABILITY + ) + profile = PoissonProfile( + mean_change=self.basis[0], + alpha=self.basis[1], + offset=self.basis[2], + tail_probability=tail_probability, + ) + support_low, support_high = profile.support() + layer_numbers = np.arange(support_low, support_high + 1) + material_occupancy = profile.occupancy(layer_numbers) + sharp_film_occupancy = (layer_numbers < 0).astype(np.float64) + return _PoissonCandidates( + profile=profile, + layer_numbers=layer_numbers, + material_occupancy=material_occupancy, + surface_occupancy=profile.surface_occupancy(layer_numbers), + # Layer n is the top filled layer when the height change is n + 1. + exposed_probability=profile.probability(layer_numbers + 1), + film_correction_occupancy=material_occupancy - sharp_film_occupancy, + tail_probability=tail_probability, + ) + + @staticmethod + def _coherent_retention(candidates): + """Return the coherent assembly's retention mask. + + Keeps every layer carrying either exposed surface or a Film + correction, and never renormalizes. This is the historical policy and + is deliberately not the height-state ensemble's. + + :param _PoissonCandidates candidates: + Unmasked profile arrays. + :returns: + Boolean mask over ``candidates.layer_numbers``. + :rtype: numpy.ndarray + :raises ValueError: + If the profile represents no material or no exposed surface. + """ + tail_probability = candidates.tail_probability + if not np.any(candidates.material_occupancy > tail_probability): + raise ValueError("Poisson surface profile has no represented material") + represented = (candidates.surface_occupancy > tail_probability) | ( + np.abs(candidates.film_correction_occupancy) > tail_probability + ) + exposed = candidates.surface_occupancy[represented] + if not np.any(exposed > tail_probability): + raise ValueError("Poisson surface profile has no represented surface") + return represented + + @staticmethod + def _height_state_selection(candidates, exact_layer_count): + """Return the retained height-state interval and its masses. + + Retains every candidate state when there are no more than + ``exact_layer_count`` of them. For a wider support it keeps a + contiguous interval containing the mode of the calculated masses, + holding at least ``exact_layer_count`` states, and expands it toward + the larger adjacent mass until the excluded mass meets the profile's + cumulative tail target. Measured CTR values never enter this choice. + + The reported masses come from the profile's own cumulative + distribution rather than from summing the candidate array, so the + retained and two excluded masses add to one however the candidate + support was truncated. + + :param _PoissonCandidates candidates: + Unmasked profile arrays. + :param int exact_layer_count: + Positive state count below which every state is retained. + :returns: + Low and high indices into ``candidates.layer_numbers``, the + normalized masses, and the raw retained, lower, and upper masses. + :rtype: tuple + :raises ValueError: + If ``exact_layer_count`` is not a positive integer or the profile + populates no height state. + """ + if int(exact_layer_count) != exact_layer_count or exact_layer_count < 1: + raise ValueError("exact_layer_count must be a positive integer") + exact_layer_count = int(exact_layer_count) + + masses = candidates.exposed_probability + populated = np.flatnonzero(masses > 0.0) + if populated.size == 0: + raise ValueError("Poisson surface profile has no populated height state") + + if populated.size <= exact_layer_count: + low, high = int(populated[0]), int(populated[-1]) + else: + # Ties resolve to the lower index. The maximum of the convolved + # masses is not floor(rate): for alpha < 1 the distribution is a + # two-point step mixture convolved with the Poisson. + low = high = int(np.argmax(masses)) + first, last = int(populated[0]), int(populated[-1]) + target = candidates.tail_probability + while low > first or high < last: + wide_enough = (high - low + 1) >= exact_layer_count + excluded = 1.0 - _retained_mass(candidates, low, high) + if wide_enough and excluded <= target: + break + below = masses[low - 1] if low > first else -1.0 + above = masses[high + 1] if high < last else -1.0 + if above >= below: + high += 1 + else: + low -= 1 + + raw_retained = _retained_mass(candidates, low, high) + if raw_retained <= 0.0: + raise ValueError("Poisson surface profile has no retained height mass") + layers = candidates.layer_numbers + excluded_lower = float( + 1.0 - candidates.profile.occupancy([layers[low]])[0] + ) + excluded_upper = float( + candidates.profile.occupancy([layers[high] + 1])[0] + ) + probabilities = masses[low : high + 1] / raw_retained + return low, high, probabilities, raw_retained, excluded_lower, excluded_upper + + def _layer_domain_geometry(self, layer_number, mat_0=None): + """Return the placement of one structural layer's three slabs. + + Shared by :meth:`createLayers` and :meth:`flat_domain_corrections` so + that layer numbering, strain, and terrace placement are derived in one + place. The matrices do not yet carry the ``below_H`` translation which + ``_translate_domains`` applies to the stored domains. + + :param int layer_number: + Structural layer offset relative to the sharp Film boundary. + :param numpy.ndarray mat_0: + Optional identity template, to avoid rebuilding it per layer. + :returns: + The three cells and their domain matrices. + :rtype: _LayerGeometry + """ + if mat_0 is None: + mat_0 = np.vstack((np.identity(3).T, np.array([0, 0, 0]))).T + n_layers_in_uc = len(self._layer_ids) + order_index = layer_number % n_layers_in_uc + cycle_index = layer_number // n_layers_in_uc + surface_uc = self.layer_ucs[order_index] + film_uc = self.film_layer_ucs[order_index] + layer_id = self.layer_order[order_index] + reference_uc = self._film_termination_ucs[float(layer_id)] + + relative_layer_position = ( + cycle_index + self.layerpos[order_index] - self.layerpos[0] + ) + layer_offset = ( + relative_layer_position + - self.underlying_film.unitcell.layerpos[layer_id] + ) + film_strain = self.underlying_film.unitcell.coherentDomainMatrix[0][2, 2] + + film_matrix = np.copy(mat_0) + film_matrix[2, 2] = film_strain + film_matrix[2, 3] = layer_offset * film_strain + + terrace_height = ( + relative_layer_position * self.underlying_film.unitcell.a[2] + ) + + surface_matrix = np.copy(mat_0) + surface_strain = self._termination_domain_strain[float(layer_id)] + surface_origin = surface_uc.layerpos[float(layer_id)] + surface_matrix[2, 2] = surface_strain + surface_matrix[2, 3] = ( + terrace_height / surface_uc.a[2] - surface_strain * surface_origin + ) + + reference_matrix = np.copy(mat_0) + reference_origin = reference_uc.layerpos[float(layer_id)] + reference_matrix[2, 2] = film_strain + reference_matrix[2, 3] = ( + terrace_height / reference_uc.a[2] - film_strain * reference_origin + ) + return _LayerGeometry( + layer_number=int(layer_number), + layer_id=float(layer_id), + surface_uc=surface_uc, + surface_matrix=surface_matrix, + film_uc=film_uc, + film_matrix=film_matrix, + reference_uc=reference_uc, + reference_matrix=reference_matrix, + ) + + def _with_below_translation(self, matrix, uc): + """Return the matrix carrying the same shift ``_translate_domains`` adds.""" + translated = np.copy(matrix) + translated[2, 3] += self.below_H / uc.a[2] + return translated + + @staticmethod + def _amplitude_for_domain(uc, matrix, occupancy, h, k, l): # noqa: E741 + """Evaluate one cell at one domain without disturbing stored domains. + + The persistent ``coherentDomainMatrix``/``coherentDomainOccupancy`` + lists are restored even when evaluation raises, so a flat-state + calculation never perturbs the coherent assembly. + """ + saved_matrix = uc.coherentDomainMatrix + saved_occupancy = uc.coherentDomainOccupancy + try: + uc.coherentDomainMatrix = [matrix] + uc.coherentDomainOccupancy = [occupancy] + return uc.F_uc(h, k, l) + finally: + uc.coherentDomainMatrix = saved_matrix + uc.coherentDomainOccupancy = saved_occupancy + + def _iter_flat_state_amplitudes(self, layer_numbers, h, k, l): # noqa: E741 + """Yield ``(layer_number, correction)`` outward from the boundary. + + The cumulative Film correction is anchored at the sharp boundary, + where it is zero, and grown upward by adding layers and downward by + removing them. Each structural layer between the boundary and the + retained interval is evaluated exactly once, and only one amplitude + array is in flight. + """ + film_occupancy = self.underlying_film.unitcell.coherentDomainOccupancy[0] + mat_0 = np.vstack((np.identity(3).T, np.array([0, 0, 0]))).T + lowest = int(layer_numbers[0]) + highest = int(layer_numbers[-1]) + prefix = np.zeros_like(np.asarray(l, dtype=np.float64), dtype=np.complex128) + + running = prefix + for layer_number in range(0, highest + 1): + geometry = self._layer_domain_geometry(layer_number, mat_0) + running = running + self._film_layer_amplitude( + geometry, film_occupancy, h, k, l + ) + if layer_number >= lowest: + yield layer_number, running + self._termination_correction( + geometry, film_occupancy, h, k, l + ) + + running = prefix + for layer_number in range(-1, lowest - 1, -1): + geometry = self._layer_domain_geometry(layer_number, mat_0) + # Layers between the boundary and the retained interval still + # have to be walked through to build the correction, but they are + # not states of the ensemble and must not be yielded. + if layer_number <= highest: + yield layer_number, running + self._termination_correction( + geometry, film_occupancy, h, k, l + ) + if layer_number > lowest: + running = running - self._film_layer_amplitude( + geometry, film_occupancy, h, k, l + ) + + def _film_layer_amplitude(self, geometry, film_occupancy, h, k, l): # noqa: E741 + """Return one Film layer slab at occupancy one.""" + return self._amplitude_for_domain( + geometry.film_uc, + self._with_below_translation(geometry.film_matrix, geometry.film_uc), + film_occupancy, + h, + k, + l, + ) + + def _termination_correction(self, geometry, film_occupancy, h, k, l): # noqa: E741 + """Return the exposed termination minus the Film material it replaces.""" + exposed = self._amplitude_for_domain( + geometry.surface_uc, + self._with_below_translation( + geometry.surface_matrix, geometry.surface_uc + ), + self._termination_domain_occupancy[geometry.layer_id], + h, + k, + l, + ) + replaced = self._amplitude_for_domain( + geometry.reference_uc, + self._with_below_translation( + geometry.reference_matrix, geometry.reference_uc + ), + -film_occupancy, + h, + k, + l, + ) + return exposed + replaced + + def flat_domain_corrections(self, h, k, l, *, exact_layer_count=10): # noqa: E741 + """Return the retained flat-height states and their corrections. + + Each correction is the occupancy-one component amplitude for one flat + surface height, relative to the same sharp Film boundary that + :meth:`F_uc` corrects. The state probability is not folded into it. + + :param numpy.ndarray h: + Reference-frame reciprocal coordinate in r.l.u. + :param numpy.ndarray k: + Reference-frame reciprocal coordinate in r.l.u. + :param numpy.ndarray l: + Reference-frame reciprocal coordinate in r.l.u. + :param int exact_layer_count: + Retain every state when the support holds no more than this many. + A policy switch, not a cap: a wider support may still retain more. + :returns: + Retained states, their normalized masses, and the excluded tails. + :rtype: FlatHeightCorrections + :raises ValueError: + If the surface is not stacked on a Film, or the profile populates + no height state. + """ + if self.underlying_film is None: + raise ValueError( + "PoissonSurface must be stacked immediately above a Film " + "before flat-height corrections can be calculated" + ) + # Pick up parameter changes through the same path as `F_uc`. + if np.any(self._basis_created != self.basis): + self.createLayers() + if ctr_accel_enabled(): + h, k, l = _ensure_contiguous( # noqa: E741 + h, k, l, testOnly=False, astype=np.float64 + ) + candidates = self._profile_candidates() + ( + low, + high, + probabilities, + raw_retained, + excluded_lower, + excluded_upper, + ) = self._height_state_selection(candidates, exact_layer_count) + layer_numbers = np.array( + candidates.layer_numbers[low : high + 1], dtype=np.int64 + ) + return FlatHeightCorrections( + layer_numbers=layer_numbers, + probabilities=np.array(probabilities, dtype=np.float64), + raw_retained_probability=float(raw_retained), + excluded_lower_probability=excluded_lower, + excluded_upper_probability=excluded_upper, + amplitude_source=lambda: self._iter_flat_state_amplitudes( + layer_numbers, h, k, l + ), + ) + def createLayers(self): """Create co-located surface and covered-Film layer domains.""" if self.underlying_film is None: @@ -2269,37 +2804,15 @@ def createLayers(self): ref.basis[:, 3] = film_basis[:, 3] + offsets n_layers_in_uc = len(self._layer_ids) - tail_probability = ( - self.profile.tail_probability - if self.profile is not None - else DEFAULT_TAIL_PROBABILITY - ) - profile = PoissonProfile( - mean_change=self.basis[0], - alpha=self.basis[1], - offset=self.basis[2], - tail_probability=tail_probability, - ) - support_low, support_high = profile.support() - layer_numbers = np.arange(support_low, support_high + 1) - material_occupancy = profile.occupancy(layer_numbers) - represented_material = np.flatnonzero( - material_occupancy > tail_probability - ) - if represented_material.size == 0: - raise ValueError("Poisson surface profile has no represented material") - surface_occupancy = profile.surface_occupancy(layer_numbers) - sharp_film_occupancy = (layer_numbers < 0).astype(np.float64) - film_correction_occupancy = material_occupancy - sharp_film_occupancy - represented = (surface_occupancy > tail_probability) | ( - np.abs(film_correction_occupancy) > tail_probability - ) - layer_numbers = layer_numbers[represented] - surface_occupancy = surface_occupancy[represented] - film_correction_occupancy = film_correction_occupancy[represented] + candidates = self._profile_candidates() + tail_probability = candidates.tail_probability + represented = self._coherent_retention(candidates) + layer_numbers = candidates.layer_numbers[represented] + surface_occupancy = candidates.surface_occupancy[represented] + film_correction_occupancy = candidates.film_correction_occupancy[ + represented + ] exposed_layers = layer_numbers[surface_occupancy > tail_probability] - if exposed_layers.size == 0: - raise ValueError("Poisson surface profile has no represented surface") top_surface_layer = exposed_layers[-1] layers_to_create = len(layer_numbers) for uc in self.layer_ucs: @@ -2318,57 +2831,25 @@ def createLayers(self): ) for layer_index, layer_number in enumerate(layer_numbers): - order_index = layer_number % n_layers_in_uc - cycle_index = layer_number // n_layers_in_uc - uc = self.layer_ucs[order_index] - film_uc = self.film_layer_ucs[order_index] - reference_uc = self._film_termination_ucs[ - float(self.layer_order[order_index]) - ] - mat_i = np.copy(mat_0) - layer_id = self.layer_order[order_index] - relative_layer_position = ( - cycle_index + self.layerpos[order_index] - self.layerpos[0] - ) - layer_offset = ( - relative_layer_position - - self.underlying_film.unitcell.layerpos[layer_id] - ) + geometry = self._layer_domain_geometry(layer_number, mat_0) - film_strain = self.underlying_film.unitcell.coherentDomainMatrix[0][2, 2] - mat_i[2, 2] = film_strain - mat_i[2, 3] = layer_offset * film_strain - - film_uc.coherentDomainMatrix.append(np.copy(mat_i)) - film_uc.coherentDomainOccupancy.append( + geometry.film_uc.coherentDomainMatrix.append(geometry.film_matrix) + geometry.film_uc.coherentDomainOccupancy.append( film_domain_occupancy * film_correction_occupancy[layer_index] ) - terrace_height = ( - relative_layer_position * self.underlying_film.unitcell.a[2] - ) - surface_matrix = np.copy(mat_0) - surface_strain = self._termination_domain_strain[float(layer_id)] - surface_origin = uc.layerpos[float(layer_id)] - surface_matrix[2, 2] = surface_strain - surface_matrix[2, 3] = ( - terrace_height / uc.a[2] - surface_strain * surface_origin + geometry.surface_uc.coherentDomainMatrix.append( + geometry.surface_matrix ) - uc.coherentDomainMatrix.append(surface_matrix) - uc.coherentDomainOccupancy.append( - self._termination_domain_occupancy[float(layer_id)] + geometry.surface_uc.coherentDomainOccupancy.append( + self._termination_domain_occupancy[geometry.layer_id] * surface_occupancy[layer_index] ) - reference_matrix = np.copy(mat_0) - reference_origin = reference_uc.layerpos[float(layer_id)] - reference_matrix[2, 2] = film_strain - reference_matrix[2, 3] = ( - terrace_height / reference_uc.a[2] - - film_strain * reference_origin + geometry.reference_uc.coherentDomainMatrix.append( + geometry.reference_matrix ) - reference_uc.coherentDomainMatrix.append(reference_matrix) - reference_uc.coherentDomainOccupancy.append( + geometry.reference_uc.coherentDomainOccupancy.append( -film_domain_occupancy * surface_occupancy[layer_index] ) diff --git a/orgui/datautils/xrayutils/test/test_CTRfilm.py b/orgui/datautils/xrayutils/test/test_CTRfilm.py new file mode 100644 index 0000000..f5025b1 --- /dev/null +++ b/orgui/datautils/xrayutils/test/test_CTRfilm.py @@ -0,0 +1,354 @@ +"""Regression tests for the Poisson surface flat-height decomposition. + +The reference for every state correction is an independently stacked ``Film`` +from ``_poisson_oracle``, which never touches ``PoissonSurface``. These cover +increment 2 of ``doc/design/incoherent_ctr_models.md``. +""" + +import copy +import unittest + +import numpy as np + +from .. import CTRfilm +from ..CTRdistributions import PoissonProfile +from . import _poisson_oracle + + +class PoissonFlatHeightMixin(unittest.TestCase): + """Shared coordinates and builders for flat-height tests.""" + + H = np.zeros(5) + K = np.zeros(5) + L = np.array([0.35, 0.8, 1.3, 1.85, 2.4]) + + def bound_surface(self, profile, n_layers=2, **keyargs): + """Return a stacked crystal, its surface, and the base width.""" + w_base = _poisson_oracle.minimum_base_width(profile) + crystal, surface = _poisson_oracle.poisson_crystal( + profile, w_base=w_base, n_layers=n_layers, **keyargs + ) + crystal.apply_stacking() + return crystal, surface, w_base + + +class TestFlatHeightCorrections(PoissonFlatHeightMixin): + """State corrections reproduce independently built flat-height crystals.""" + + def test_state_corrections_match_the_flat_height_oracle(self): + """Adding one state correction to the base crystal gives that height. + + The base crystal is the same sharp Film boundary the correction is + defined against, so ``common + correction`` must reproduce a crystal + independently stacked at that flat height. + """ + cases = { + "integer growth": PoissonProfile(2.0, alpha=0.0), + "fractional step": PoissonProfile(1.5, alpha=0.0), + "negative offset": PoissonProfile(0.0, alpha=0.0, offset=-2.25), + "poisson growth": PoissonProfile(1.5, alpha=0.6), + "poisson etching": PoissonProfile(-1.5, alpha=0.6), + "step and poisson": PoissonProfile(2.3, alpha=0.4, offset=0.4), + } + for label, profile in cases.items(): + with self.subTest(case=label): + profile.tail_probability = 1e-14 + _, surface, w_base = self.bound_surface(profile) + common = _poisson_oracle.flat_height_crystal( + -1, w_base=w_base + ).F(self.H, self.K, self.L) + states = surface.flat_domain_corrections( + self.H, self.K, self.L + ) + self.assertGreater(states.layer_numbers.size, 0) + for state in states.iter_states(): + oracle = _poisson_oracle.flat_height_crystal( + state.layer_number, w_base=w_base + ).F(self.H, self.K, self.L) + np.testing.assert_allclose( + common + state.amplitude, oracle, rtol=1e-12 + ) + + def test_multi_layer_cycles_and_termination_banks(self): + """Layer cycles and an explicit termination bank keep the agreement.""" + profile = PoissonProfile(1.5, alpha=0.5, tail_probability=1e-14) + for n_layers in (1, 2, 3): + with self.subTest(n_layers=n_layers): + _, surface, w_base = self.bound_surface( + profile, n_layers=n_layers + ) + common = _poisson_oracle.flat_height_crystal( + -1, w_base=w_base, n_layers=n_layers + ).F(self.H, self.K, self.L) + states = surface.flat_domain_corrections( + self.H, self.K, self.L + ) + for state in states.iter_states(): + oracle = _poisson_oracle.flat_height_crystal( + state.layer_number, + w_base=w_base, + n_layers=n_layers, + ).F(self.H, self.K, self.L) + np.testing.assert_allclose( + common + state.amplitude, oracle, rtol=1e-12 + ) + + def test_copied_surfaces_still_produce_corrections(self): + """``copy.deepcopy`` is performed once per fit and must survive it.""" + profile = PoissonProfile(1.5, alpha=0.5, tail_probability=1e-14) + crystal, surface, w_base = self.bound_surface(profile) + original = surface.flat_domain_corrections( + self.H, self.K, self.L + ).as_array() + + copied_crystal = copy.deepcopy(crystal) + copied_crystal.apply_stacking() + copied = copied_crystal.uc_surface_list[1].flat_domain_corrections( + self.H, self.K, self.L + ).as_array() + np.testing.assert_allclose(copied, original, rtol=1e-12) + + def test_parameter_changes_are_picked_up(self): + """Corrections follow the basis through the ``F_uc`` sync path.""" + profile = PoissonProfile(1.5, alpha=0.0, tail_probability=1e-14) + _, surface, _ = self.bound_surface(profile) + before = surface.flat_domain_corrections(self.H, self.K, self.L) + + surface.basis[0] = 3.0 + after = surface.flat_domain_corrections(self.H, self.K, self.L) + self.assertNotEqual( + before.layer_numbers.tolist(), after.layer_numbers.tolist() + ) + + +class TestFlatHeightSelection(PoissonFlatHeightMixin): + """The ensemble selection policy and its reported masses.""" + + def test_interior_bins_agree_with_the_shifted_probability(self): + """``surface_occupancy(n)`` is ``probability(n + 1)`` except at the top. + + Layer ``n`` is the top filled layer exactly when the height change is + ``n + 1``. The two expressions therefore agree on every interior bin, + and deliberately disagree at the terminal one, where + ``surface_occupancy`` folds in the whole upper tail. Both halves are + asserted so a refactor cannot silently exchange them. + """ + profile = PoissonProfile(1.5, alpha=0.6, tail_probability=1e-6) + _, surface, _ = self.bound_surface(profile) + candidates = surface._profile_candidates() + + np.testing.assert_allclose( + candidates.surface_occupancy[:-1], + candidates.exposed_probability[:-1], + atol=1e-15, + ) + self.assertGreater( + abs( + candidates.surface_occupancy[-1] + - candidates.exposed_probability[-1] + ), + 0.0, + ) + + def test_the_two_retention_policies_select_different_sets(self): + """Coherent retention keeps Film-correction layers the ensemble drops. + + A single deterministic grown state exposes only its own layer, but the + coherent assembly must also keep the Film layer it added underneath. + Merging the two policies into one mask would silently give one path + the other's cutoff. + """ + profile = PoissonProfile(2.0, alpha=0.0, tail_probability=1e-14) + _, surface, _ = self.bound_surface(profile) + candidates = surface._profile_candidates() + + coherent = candidates.layer_numbers[ + surface._coherent_retention(candidates) + ] + low, high, _, _, _, _ = surface._height_state_selection(candidates, 10) + ensemble = candidates.layer_numbers[low : high + 1] + + self.assertEqual(ensemble.tolist(), [1]) + self.assertEqual(coherent.tolist(), [0, 1]) + + def test_exact_layer_count_is_a_policy_not_a_cap(self): + """Below the count every state is kept; above it the mass decides.""" + narrow = PoissonProfile(1.5, alpha=0.0, tail_probability=1e-14) + _, surface, _ = self.bound_surface(narrow) + states = surface.flat_domain_corrections( + self.H, self.K, self.L, exact_layer_count=10 + ) + self.assertEqual(states.layer_numbers.size, 2) + + wide = PoissonProfile(6.0, alpha=1.0, tail_probability=1e-6) + _, surface, _ = self.bound_surface(wide) + candidates = surface._profile_candidates() + populated = int(np.sum(candidates.exposed_probability > 0.0)) + self.assertGreater(populated, 10) + + # More than ten states are retained because the probability-mass + # criterion demands them, and asking for fewer does not cap it. + for requested in (3, 10): + with self.subTest(exact_layer_count=requested): + states = surface.flat_domain_corrections( + self.H, self.K, self.L, exact_layer_count=requested + ) + self.assertGreater(states.layer_numbers.size, 10) + self.assertLessEqual( + states.excluded_probability, wide.tail_probability + ) + + everything = surface.flat_domain_corrections( + self.H, self.K, self.L, exact_layer_count=populated + 5 + ) + self.assertEqual(everything.layer_numbers.size, populated) + + def test_masses_are_normalized_once_and_add_to_one(self): + """Retained and excluded masses partition the whole distribution.""" + for tail in (1e-14, 1e-3): + with self.subTest(tail_probability=tail): + profile = PoissonProfile(2.0, alpha=1.0, tail_probability=tail) + _, surface, _ = self.bound_surface(profile) + states = surface.flat_domain_corrections( + self.H, self.K, self.L + ) + self.assertAlmostEqual( + float(states.probabilities.sum()), 1.0, places=12 + ) + self.assertAlmostEqual( + states.raw_retained_probability + + states.excluded_lower_probability + + states.excluded_upper_probability, + 1.0, + places=12, + ) + self.assertAlmostEqual( + states.excluded_probability, + states.excluded_lower_probability + + states.excluded_upper_probability, + places=15, + ) + + def test_excluded_mass_is_not_an_amplitude_error_bound(self): + """A looser tail costs more amplitude error than dimensionless mass. + + This is why the design record refuses to use the excluded probability + directly as a reconstruction tolerance. + """ + profile = PoissonProfile(2.0, alpha=1.0, tail_probability=1e-3) + _, surface, _ = self.bound_surface(profile) + direct = surface.F_uc(self.H, self.K, self.L) + states = surface.flat_domain_corrections(self.H, self.K, self.L) + rebuilt = sum( + state.probability * state.amplitude + for state in states.iter_states() + ) + relative = np.max(np.abs(direct - rebuilt)) / np.max(np.abs(direct)) + self.assertGreater(states.excluded_probability, 0.0) + self.assertGreater(relative, states.excluded_probability) + + def test_tight_tail_reconstructs_the_coherent_amplitude(self): + """With a negligible tail the state sum reproduces ``F_uc``.""" + for label, profile in { + "growth": PoissonProfile(1.5, alpha=0.6, tail_probability=1e-14), + "etching": PoissonProfile(-1.5, alpha=0.6, tail_probability=1e-14), + }.items(): + with self.subTest(case=label): + _, surface, _ = self.bound_surface(profile) + direct = surface.F_uc(self.H, self.K, self.L) + states = surface.flat_domain_corrections( + self.H, self.K, self.L + ) + rebuilt = sum( + state.probability * state.amplitude + for state in states.iter_states() + ) + np.testing.assert_allclose(rebuilt, direct, rtol=1e-11) + + +class TestFlatHeightContract(PoissonFlatHeightMixin): + """API behavior of the returned result.""" + + def test_evaluation_leaves_the_coherent_assembly_untouched(self): + """Flat states must not perturb the stored coherent domains.""" + profile = PoissonProfile(1.5, alpha=0.5, tail_probability=1e-14) + _, surface, _ = self.bound_surface(profile) + before = surface.F_uc(self.H, self.K, self.L) + snapshot = [ + ( + [np.copy(matrix) for matrix in uc.coherentDomainMatrix], + list(uc.coherentDomainOccupancy), + ) + for uc in surface.layer_ucs + surface.film_layer_ucs + ] + + surface.flat_domain_corrections( + self.H, self.K, self.L + ).as_array() + + for uc, (matrices, occupancies) in zip( + surface.layer_ucs + surface.film_layer_ucs, snapshot + ): + self.assertEqual(len(uc.coherentDomainMatrix), len(matrices)) + for stored, expected in zip(uc.coherentDomainMatrix, matrices): + np.testing.assert_array_equal(stored, expected) + np.testing.assert_allclose( + uc.coherentDomainOccupancy, occupancies + ) + np.testing.assert_array_equal( + surface.F_uc(self.H, self.K, self.L), before + ) + + def test_as_array_matches_the_streamed_states(self): + """The diagnostic materialization agrees with the streaming path.""" + profile = PoissonProfile(1.5, alpha=0.5, tail_probability=1e-14) + _, surface, _ = self.bound_surface(profile) + states = surface.flat_domain_corrections(self.H, self.K, self.L) + + streamed = { + state.layer_number: state.amplitude + for state in states.iter_states() + } + materialized = states.as_array() + self.assertEqual( + materialized.shape, (states.layer_numbers.size, self.L.size) + ) + for index, layer in enumerate(states.layer_numbers): + np.testing.assert_array_equal( + materialized[index], streamed[int(layer)] + ) + + def test_reported_arrays_are_read_only(self): + """Callers must not be able to edit the reported selection.""" + profile = PoissonProfile(1.5, alpha=0.5, tail_probability=1e-14) + _, surface, _ = self.bound_surface(profile) + states = surface.flat_domain_corrections(self.H, self.K, self.L) + with self.assertRaises(ValueError): + states.layer_numbers[0] = 0 + with self.assertRaises(ValueError): + states.probabilities[0] = 0.0 + + def test_requires_a_bound_film(self): + """A surface which is not stacked has no boundary to correct.""" + surface = CTRfilm.PoissonSurface( + _poisson_oracle.layered_cell(2, "film"), + profile=PoissonProfile(1.5, alpha=0.5), + name="surface", + ) + with self.assertRaisesRegex(ValueError, "stacked immediately above"): + surface.flat_domain_corrections(self.H, self.K, self.L) + + def test_rejects_an_invalid_state_count(self): + """``exact_layer_count`` is a positive integer policy setting.""" + profile = PoissonProfile(1.5, alpha=0.5, tail_probability=1e-14) + _, surface, _ = self.bound_surface(profile) + for bad in (0, -1, 2.5): + with self.subTest(exact_layer_count=bad): + with self.assertRaisesRegex(ValueError, "positive integer"): + surface.flat_domain_corrections( + self.H, self.K, self.L, exact_layer_count=bad + ) + + +if __name__ == "__main__": + unittest.main() From 95a15d8834854cd7847fa1d168e66ff6cf286efa Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Fri, 11 Sep 2026 18:50:09 -0400 Subject: [PATCH 09/13] feat(phys): add the incoherent CTR model contract Wrappers present the existing fit API with their own parameters ahead of the wrapped crystal's, and mix complete coherent state amplitudes. The registry lets a saved type name resolve without importing a class path. --- orgui/datautils/xrayutils/CTRcalc.py | 23 +- orgui/datautils/xrayutils/CTRincoherent.py | 836 ++++++++++++++++++ .../xrayutils/test/test_CTRincoherent.py | 667 ++++++++++++++ 3 files changed, 1522 insertions(+), 4 deletions(-) create mode 100644 orgui/datautils/xrayutils/CTRincoherent.py create mode 100644 orgui/datautils/xrayutils/test/test_CTRincoherent.py diff --git a/orgui/datautils/xrayutils/CTRcalc.py b/orgui/datautils/xrayutils/CTRcalc.py index 903819c..25ae1a3 100644 --- a/orgui/datautils/xrayutils/CTRcalc.py +++ b/orgui/datautils/xrayutils/CTRcalc.py @@ -450,13 +450,28 @@ def _iter_scaled_contributions(self, index, hkl, evaluator): electrons per the component's own lateral cell. """ uc = self.uc_surface_list[index] + for matrix, scale in self._component_domain_factors(index): + hkl_n = np.dot(matrix, hkl) + yield scale * evaluator(uc, hkl_n[0], hkl_n[1], hkl_n[2]) + + def _component_domain_factors(self, index): + """Yield the transform and combined scale of each outer domain. + + The scale is ``reference_area / component_area * occupancy * weight``. + This is the single definition of that product: an incoherent model + combining its own per-state amplitudes must use it rather than + rebuild it, or the two paths can drift apart. + + :param int index: + Position in :attr:`uc_surface_list`. + :returns: + Generator of ``(matrix, scale)`` pairs in domain order. + """ + uc = self.uc_surface_list[index] area_scale = self.reference_area / uc.uc_area weight = self.weights[index] for matrix, occup in self.domains[index]: - hkl_n = np.dot(matrix, hkl) - yield area_scale * occup * weight * evaluator( - uc, hkl_n[0], hkl_n[1], hkl_n[2] - ) + yield matrix, area_scale * occup * weight def evaluate_kinematic(self, harray, karray, Larray): """Return the decomposed kinematical evaluation behind :meth:`F`. diff --git a/orgui/datautils/xrayutils/CTRincoherent.py b/orgui/datautils/xrayutils/CTRincoherent.py new file mode 100644 index 0000000..c393368 --- /dev/null +++ b/orgui/datautils/xrayutils/CTRincoherent.py @@ -0,0 +1,836 @@ +# /*########################################################################## +# +# Copyright (c) 2020-2025 Timo Fuchs +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ +"""Incoherent kinematical CTR models. + +An incoherent model wraps one coherent :class:`CTRcalc.SXRDCrystal` and +returns a squared structure factor for a mixed state. Constructing a wrapper +is the only way to opt in: ``SXRDCrystal.F`` and every existing calculation +stay fully coherent. + +Quantities, following ``doc/design/incoherent_ctr_models.md``: + +- ``F`` is the coherent complex amplitude and exists only on the crystal. A + mixed state has no unique complex amplitude, so no wrapper defines ``F``. +- ``F2`` is ``abs(F) ** 2`` in the squared units of ``F``, nominally electrons + squared per reference lateral cell. It is not detector counts and not a + reflectivity. + +The incoherent average acts on the amplitude of the *complete* coherent +crystal state, not on the isolated surface correction, so bulk-surface and +Film-surface interference stays inside each domain. +""" + +__author__ = "Timo Fuchs" +__copyright__ = "Copyright 2020-2025 Timo Fuchs" +__license__ = "MIT License" +__maintainer__ = "Timo Fuchs" +__email__ = "tfuchs@cornell.edu" + +import inspect +from abc import ABC, abstractmethod +from dataclasses import dataclass +from types import MappingProxyType + +import numpy as np + +from .CTRcalc import SXRDCrystal +from .CTRutil import LinearFitFunctions + +__all__ = [ + "CoherentStateEnsembleModel", + "IncoherentF2Model", + "IncoherentModel", + "IncoherentModelInfo", + "IncoherentState", + "KinematicIncoherentContext", + "available_incoherent_models", + "create_incoherent_model", + "create_incoherent_model_from_config", + "register_incoherent_model", +] + + +@dataclass(frozen=True) +class IncoherentState: + """One complete coherent state of the wrapped crystal. + + :param int index: + Position of the state in its source's own ordering. + :param float probability: + Normalized weight of this state in the mixture. + :param amplitude: + Complete crystal amplitude when this state covers the coherent patch, + in electrons per reference lateral cell. + """ + + index: int + probability: float + amplitude: object + + +class KinematicIncoherentContext: + """Lazy, evaluation-local decomposition for one HKL request. + + The coherent decomposition is calculated at most once per context, so the + semi-infinite bulk and every common component are evaluated once no matter + how many states a model streams. Nothing is cached beyond the context, + which is what keeps the cache free of invalidation hazards: optimizer + parameters, component weights, domain transforms, stacking, attenuation, + and the coordinates themselves all change between requests. + + :param CTRcalc.SXRDCrystal crystal: + The wrapped coherent crystal. + :param numpy.ndarray h: + Reference-frame reciprocal coordinate in r.l.u. + :param numpy.ndarray k: + Reference-frame reciprocal coordinate in r.l.u. + :param numpy.ndarray l: + Reference-frame reciprocal coordinate in r.l.u. + """ + + def __init__(self, crystal, h, k, l): # noqa: E741 + self._crystal = crystal + self._h = h + self._k = k + self._l = l + self._coherent = None + self._hkl = None + + @property + def crystal(self): + """Return the wrapped coherent crystal.""" + return self._crystal + + @property + def coordinates(self): + """Return the requested ``(h, k, l)`` coordinates in r.l.u.""" + return self._h, self._k, self._l + + @property + def coherent(self): + """Return the coherent decomposition, evaluating it on first use. + + :rtype: CTRcalc.KinematicAmplitudeResult + """ + if self._coherent is None: + self._coherent = self._crystal.evaluate_kinematic( + self._h, self._k, self._l + ) + return self._coherent + + def component(self, name): + """Return one component's amplitude record, selected by stable name. + + :param str name: + Component name on the wrapped crystal. + :rtype: CTRcalc.KinematicComponentAmplitude + """ + return self.coherent.component(name) + + def common_amplitude(self, name): + """Return the crystal amplitude with one component removed. + + This is the part shared by every state of that component: the bulk, + the underlying Film, and all other components. + + :param str name: + Component name on the wrapped crystal. + :returns: + Complex amplitude in electrons per reference lateral cell. + """ + return self.coherent.total - self.component(name).amplitude + + def _stacked_hkl(self): + """Return the stacked ``(3, N)`` coordinates used by domain transforms.""" + if self._hkl is None: + self._hkl = np.vstack((self._h, self._k, self._l)) + return self._hkl + + def iter_component_states(self, name, state_evaluator): + """Yield one complete state amplitude of a component at a time. + + ``state_evaluator`` is called once per outer coherent domain, at that + domain's transformed coordinates, and must return an object exposing + ``iter_states()`` over items carrying ``probability`` and + ``amplitude``. Calling it at the transformed coordinates is required: + an outer domain changes where a component state has to be evaluated, + so an amplitude precomputed only at the original coordinates would be + placed wrongly. + + Every domain must report the same states in the same order. Only the + amplitudes may differ between domains, because only the coordinates + do. + + :param str name: + Component name on the wrapped crystal. + :param callable state_evaluator: + Called as ``state_evaluator(component, h, k, l)``. + :returns: + Generator of :class:`IncoherentState`. + :raises ValueError: + If the domains disagree on the state metadata. + """ + part = self.component(name) + common = self.common_amplitude(name) + component = self._crystal.uc_surface_list[part.index] + hkl = self._stacked_hkl() + + factors = list(self._crystal._component_domain_factors(part.index)) + sources = [] + for matrix, scale in factors: + transformed = np.dot(matrix, hkl) + sources.append( + ( + scale, + state_evaluator( + component, + transformed[0], + transformed[1], + transformed[2], + ), + ) + ) + + iterators = [source.iter_states() for _, source in sources] + scales = [scale for scale, _ in sources] + for index, group in enumerate(zip(*iterators)): + reference = group[0] + for other in group[1:]: + same_layer = getattr(other, "layer_number", None) == getattr( + reference, "layer_number", None + ) + if not same_layer or not np.isclose( + other.probability, reference.probability + ): + raise ValueError( + "Coherent domains of component " + f"{name!r} disagree on state {index}: only the state " + "amplitudes may depend on the domain transform" + ) + delta = None + for scale, state in zip(scales, group): + scaled = scale * state.amplitude + delta = scaled if delta is None else delta + scaled + yield IncoherentState( + index=index, + probability=float(reference.probability), + amplitude=common if delta is None else common + delta, + ) + + +class IncoherentModel(LinearFitFunctions, ABC): + """Quantity-neutral wrapper around one coherent crystal. + + The wrapper presents the fit API the optimizer already consumes, with its + own local parameters first and the wrapped crystal's parameters as the + tail:: + + [incoherent-model-local parameters] [wrapped coherent-crystal parameters] + + ``LinearFitFunctions`` already implements ``getInitialParameters``, + ``getStartParamAndLimits``, ``setFitErrors``, ``getFitErrors``, + ``fitparnames``, ``priors`` and ``parameter_list`` as **local-only** + methods. Every one of them is overridden here as a composite. Inheriting + one by accident does not raise; it silently returns a vector one block + short, which is the failure this contract exists to prevent. + + :param CTRcalc.SXRDCrystal crystal: + The coherent model to wrap. + :param str name: + Wrapper name, used as the default prefix for local parameter names. + """ + + model_type = None + supported_forward_models = frozenset({"kinematical"}) + output_quantity = None + + parameterLookup = {} + parameterLookup_inv = {} + + def __init__(self, crystal, *, name="incoherent"): + super().__init__() + if not isinstance(crystal, SXRDCrystal): + raise TypeError( + "An incoherent model wraps one SXRDCrystal, not " + f"{type(crystal).__name__}" + ) + self._coherent_model = crystal + self.name = name + + @property + def coherent_model(self): + """Return the wrapped coherent crystal.""" + return self._coherent_model + + # -- parameter block sizes ------------------------------------------ + + @property + def n_local_parameters(self): + """Return the number of wrapper-local fit parameters.""" + return len(LinearFitFunctions.parameter_list(self)) + + @property + def n_coherent_parameters(self): + """Return the number of wrapped-crystal fit parameters.""" + return len(self._coherent_model.fitparnames) + + def _split(self, values, what): + """Split a full vector into its local and coherent blocks.""" + values = np.asarray(values) + expected = self.n_local_parameters + self.n_coherent_parameters + if values.shape[0] != expected: + raise ValueError( + f"{what} has {values.shape[0]} entries but this model has " + f"{expected}: {self.n_local_parameters} local and " + f"{self.n_coherent_parameters} coherent" + ) + local = self.n_local_parameters + return values[:local], values[local:] + + # -- composite fit API ---------------------------------------------- + + def getStartParamAndLimits(self, force_recalculate=False): # noqa: N802 + """Return start values and bounds, local block first.""" + local = LinearFitFunctions.getStartParamAndLimits( + self, force_recalculate + ) + coherent = self._coherent_model.getStartParamAndLimits( + force_recalculate + ) + return tuple( + np.concatenate( + ( + np.asarray(a, dtype=np.float64), + np.asarray(b, dtype=np.float64), + ) + ) + for a, b in zip(local, coherent) + ) + + def getInitialParameters(self, force_recalculate=False): # noqa: N802 + """Return the full start vector, local block first.""" + return self.getStartParamAndLimits(force_recalculate)[0] + + def setParameters(self, values): # noqa: N802 + """Set the full parameter vector, local block first. + + The complete input is validated before either block is written, so a + wrong length leaves the model untouched rather than half updated. + + :param values: + Full parameter vector. + :raises ValueError: + If the length does not match the model. + """ + local, coherent = self._split(values, "parameter vector") + if local.size: + LinearFitFunctions.setFitParameters(self, local) + self._coherent_model.setParameters(coherent) + + def setFitParameters(self, values): # noqa: N802 + """Refuse the local-only setter inherited from ``LinearFitFunctions``. + + :raises NotImplementedError: + Always. ``setParameters`` is the composite entry point; the + inherited name would write only the local block and silently + discard the coherent tail. + """ + raise NotImplementedError( + f"{type(self).__name__}.setFitParameters would set only the " + "wrapper-local block. Use setParameters, which covers the local " + "parameters and the wrapped crystal." + ) + + def setLimits(self, limits): # noqa: N802 + """Set bounds for both blocks, local block first.""" + local, coherent = self._split(limits, "bounds array") + if local.size: + LinearFitFunctions.setLimits(self, local) + self._coherent_model.setLimits(coherent) + + def setFitErrors(self, errors): # noqa: N802 + """Set or clear errors on both blocks.""" + if errors is None: + LinearFitFunctions.setFitErrors(self, None) + self._coherent_model.setFitErrors(None) + return + local, coherent = self._split(errors, "error vector") + if local.size: + LinearFitFunctions.setFitErrors(self, local) + self._coherent_model.setFitErrors(coherent) + + def getFitErrors(self): # noqa: N802 + """Return errors for both blocks, local block first. + + A block with no fit parameters is skipped rather than queried: the + existing implementations raise when no errors are set, which for an + empty block would be a meaningless failure. A non-empty block whose + errors were never set still raises, as it does today. + """ + blocks = [] + if self.n_local_parameters: + blocks.append(np.asarray(LinearFitFunctions.getFitErrors(self))) + else: + blocks.append(np.array([], dtype=np.float64)) + if self.n_coherent_parameters: + blocks.append(np.asarray(self._coherent_model.getFitErrors())) + else: + blocks.append(np.array([], dtype=np.float64)) + return np.concatenate(blocks) + + @property + def fitparnames(self): + """Return parameter names, local block first.""" + return list(LinearFitFunctions.fitparnames.fget(self)) + list( + self._coherent_model.fitparnames + ) + + @property + def priors(self): + """Return priors, local block first.""" + return list(LinearFitFunctions.priors.fget(self)) + list( + self._coherent_model.priors + ) + + def parameter_list(self): + """Return parameter records, local block first.""" + return list(LinearFitFunctions.parameter_list(self)) + list( + self._coherent_model.parameter_list() + ) + + # -- composite serialization ---------------------------------------- + + def _local_parameters_to_dict(self): + """Return only the wrapper-local parameter subtree.""" + return LinearFitFunctions.parametersToDict(self) + + def _local_parameters_from_dict(self, data, override_values=True): + """Restore only the wrapper-local parameter subtree.""" + LinearFitFunctions.parametersFromDict(self, data, override_values) + + def parametersToDict(self): # noqa: N802 + """Return the local subtree and the wrapped crystal's subtree.""" + return { + "local": self._local_parameters_to_dict(), + "coherent": self._coherent_model.parametersToDict(), + } + + def parametersFromDict(self, data, override_values=True): # noqa: N802 + """Restore both subtrees.""" + self._local_parameters_from_dict(data["local"], override_values) + self._coherent_model.parametersFromDict( + data["coherent"], override_values + ) + + def clearParameters(self): # noqa: N802 + """Clear both blocks.""" + LinearFitFunctions.clearParameters(self) + self._coherent_model.clearParameters() + + # -- validation and configuration ----------------------------------- + + def validate(self, *, forward_model): + """Check that this model can serve the selected forward model. + + :param str forward_model: + Name of the calculation the optimizer will run. + :raises ValueError: + If the forward model is unsupported or the wrapped crystal does + not satisfy the model's own requirements. + """ + if forward_model not in self.supported_forward_models: + supported = ", ".join(sorted(self.supported_forward_models)) + raise ValueError( + f"{self.model_type or type(self).__name__} does not support " + f"the {forward_model} forward model; supported: {supported}" + ) + self._validate_model(forward_model) + + def _validate_model(self, forward_model): + """Model-specific validation hook. Override as needed.""" + + def _config_settings(self): + """Return the non-basis settings to persist. Override as needed.""" + return {} + + def to_config(self): + """Return a plain dictionary describing this model. + + The local parameter subtree is the only stored source for local basis + values, whether they are fitted or fixed; settings never duplicate + them. The wrapped crystal is not included: the factory receives it + separately, so the coherent structure keeps its own canonical file. + + :rtype: dict + """ + if not self.model_type: + raise ValueError( + f"{type(self).__name__} has no model_type and cannot be " + "serialized; register it first" + ) + return { + "type": self.model_type, + "settings": dict(self._config_settings()), + "parameters": self._local_parameters_to_dict(), + } + + +class IncoherentF2Model(IncoherentModel, ABC): + """Kinematical incoherent model returning a squared structure factor.""" + + output_quantity = "F2" + + def F2(self, h, k, l): # noqa: N802,E741 + """Return the mixed squared structure factor. + + :param numpy.ndarray h: + Reference-frame reciprocal coordinate in r.l.u. + :param numpy.ndarray k: + Reference-frame reciprocal coordinate in r.l.u. + :param numpy.ndarray l: + Reference-frame reciprocal coordinate in r.l.u. + :returns: + Real, finite, nonnegative squared structure factor in the squared + units of ``SXRDCrystal.F``. + :rtype: numpy.ndarray + :raises ValueError: + If the model produced a negative or non-finite result. + """ + context = KinematicIncoherentContext(self._coherent_model, h, k, l) + value = np.asarray(self._evaluate_F2(context), dtype=np.float64) + if not np.all(np.isfinite(value)): + raise ValueError( + f"{type(self).__name__} produced a non-finite F2" + ) + if np.any(value < 0.0): + raise ValueError( + f"{type(self).__name__} produced a negative F2" + ) + return value + + @abstractmethod + def _evaluate_F2(self, context): # noqa: N802 + """Return the mixed ``F2`` for one evaluation context.""" + + +class CoherentStateEnsembleModel(IncoherentF2Model, ABC): + """Mix the complete coherent amplitudes of a set of states. + + A subclass supplies normalized probabilities and complete state + amplitudes; this class owns the validated coherent, incoherent, and + partially coherent sums. It is deliberately generic: it neither imports + nor tests for any particular surface model. + + The endpoints are + + .. code-block:: text + + F2_coherent = abs(sum_n p_n A_n) ** 2 + F2_incoherent = sum_n p_n abs(A_n) ** 2 + + and the model returns their convex interpolation. The ``kappa = 0`` + endpoint uses the live coherent evaluation already cached in the context, + not the state-sum reconstruction, so a finite state support cannot perturb + it. + """ + + #: Absolute tolerance on the sum of the streamed probabilities. + probability_tolerance = 1e-9 + + @property + @abstractmethod + def incoherent_fraction(self): + """Return the mixing fraction ``kappa`` in ``[0, 1]``.""" + + @abstractmethod + def _iter_states(self, context): + """Yield :class:`IncoherentState` for one evaluation context.""" + + def _validated_fraction(self): + """Return the mixing fraction after range and finiteness checks.""" + kappa = float(self.incoherent_fraction) + if not np.isfinite(kappa) or not 0.0 <= kappa <= 1.0: + raise ValueError( + f"incoherent_fraction must be finite and within [0, 1], " + f"got {kappa}" + ) + return kappa + + def accumulate_states(self, context): + """Stream the states once, returning the incoherent sum and diagnostics. + + Exactly one state amplitude is held at a time: the squared sum and the + diagnostic coherent sum are accumulated in place. + + :param KinematicIncoherentContext context: + Evaluation context. + :returns: + Incoherent ``F2``, the reconstructed coherent amplitude, and the + total streamed probability. + :rtype: tuple + :raises ValueError: + If a probability is invalid or the masses do not sum to one. + """ + incoherent = None + reconstructed = None + mass = 0.0 + count = 0 + for state in self._iter_states(context): + probability = float(state.probability) + if not np.isfinite(probability) or probability < 0.0: + raise ValueError( + f"State {count} has an invalid probability {probability}" + ) + amplitude = np.asarray(state.amplitude) + if not np.all(np.isfinite(amplitude)): + raise ValueError(f"State {count} has a non-finite amplitude") + squared = probability * np.abs(amplitude) ** 2 + weighted = probability * amplitude + incoherent = squared if incoherent is None else incoherent + squared + reconstructed = ( + weighted if reconstructed is None else reconstructed + weighted + ) + mass += probability + count += 1 + if count == 0: + raise ValueError( + f"{type(self).__name__} produced no states to mix" + ) + if abs(mass - 1.0) > self.probability_tolerance: + raise ValueError( + f"State probabilities sum to {mass}, not one; they must be " + "normalized over the retained states exactly once" + ) + return incoherent, reconstructed, mass + + def _evaluate_F2(self, context): # noqa: N802 + """Interpolate between the live coherent and incoherent endpoints.""" + kappa = self._validated_fraction() + coherent = np.abs(context.coherent.total) ** 2 + if kappa == 0.0: + # Identical to the wrapped crystal's F2, without streaming any + # state or triggering a second bulk and Film evaluation. + return coherent + incoherent, _, _ = self.accumulate_states(context) + return (1.0 - kappa) * coherent + kappa * incoherent + + def coherence_reconstruction_error(self, h, k, l): # noqa: E741 + """Return the diagnostic gap between the two coherent expressions. + + The state-sum reconstruction of the coherent amplitude is a + diagnostic against the finite-support truncation; it never replaces + the live coherent endpoint. The excluded probability alone is not an + amplitude bound, so compare this against a tolerance justified by an + extended-support reference. + + :returns: + Maximum absolute difference between the live coherent amplitude + and the probability-weighted state sum. + :rtype: float + """ + context = KinematicIncoherentContext(self._coherent_model, h, k, l) + _, reconstructed, _ = self.accumulate_states(context) + return float( + np.max(np.abs(context.coherent.total - reconstructed)) + ) + + +@dataclass(frozen=True) +class IncoherentModelInfo: + """Public description of one registered incoherent model. + + :param str key: + Stable persistence key, taken from the class ``model_type``. + :param type model_class: + The registered class. + :param str description: + One-line summary from the class docstring. + :param frozenset supported_forward_models: + Forward models the class accepts. + :param str output_quantity: + Name of the quantity the class returns, such as ``"F2"``. + """ + + key: str + model_class: type + description: str + supported_forward_models: frozenset + output_quantity: str + + +_REGISTRY = {} +_BUILTIN_KEYS = set() + + +def register_incoherent_model(cls=None, *, builtin=False): + """Register one incoherent model class under its ``model_type``. + + The class attribute is the only key source, so a saved configuration and + the class cannot disagree. Usable directly or as a decorator. + + :param type cls: + Concrete :class:`IncoherentModel` subclass. + :param bool builtin: + Mark the entry as shipped with the package, which then refuses to be + replaced by a later registration. + :returns: + The class, so this works as a decorator. + :raises TypeError: + If the class is not a concrete ``IncoherentModel`` subclass. + :raises ValueError: + If the key is missing, already used, or the class has no output + quantity. + """ + if cls is None: + def decorator(inner): + return register_incoherent_model(inner, builtin=builtin) + + return decorator + + if not (isinstance(cls, type) and issubclass(cls, IncoherentModel)): + raise TypeError( + "Only IncoherentModel subclasses can be registered, not " + f"{cls!r}" + ) + if inspect.isabstract(cls): + raise TypeError( + f"{cls.__name__} is abstract and cannot be registered" + ) + key = getattr(cls, "model_type", None) + if not isinstance(key, str) or not key: + raise ValueError( + f"{cls.__name__} needs a non-empty model_type to be registered" + ) + if not getattr(cls, "output_quantity", None): + raise ValueError( + f"{cls.__name__} needs a stable output_quantity to be registered" + ) + if key in _REGISTRY: + if key in _BUILTIN_KEYS: + raise ValueError( + f"{key!r} is a built-in incoherent model and cannot be " + "replaced" + ) + raise ValueError(f"An incoherent model is already registered as {key!r}") + + summary = (inspect.getdoc(cls) or "").strip().splitlines() + _REGISTRY[key] = IncoherentModelInfo( + key=key, + model_class=cls, + description=summary[0] if summary else "", + supported_forward_models=frozenset(cls.supported_forward_models), + output_quantity=cls.output_quantity, + ) + if builtin: + _BUILTIN_KEYS.add(key) + return cls + + +def unregister_incoherent_model(key): + """Remove one non-built-in registration. + + :param str key: + Registry key. + :raises KeyError: + If the key is not registered. + :raises ValueError: + If the key names a built-in model. + """ + if key not in _REGISTRY: + raise KeyError(f"No incoherent model registered as {key!r}") + if key in _BUILTIN_KEYS: + raise ValueError(f"{key!r} is a built-in incoherent model") + del _REGISTRY[key] + + +def available_incoherent_models(): + """Return a read-only mapping from registry key to model description. + + :rtype: types.MappingProxyType + """ + return MappingProxyType(dict(_REGISTRY)) + + +def _registered(key): + """Return one registry record, or raise with the available keys.""" + try: + return _REGISTRY[key] + except (KeyError, TypeError): + available = ", ".join(sorted(_REGISTRY)) or "none" + raise ValueError( + f"Unknown incoherent model {key!r}. Registered: {available}" + ) from None + + +def create_incoherent_model(key, crystal, **keyargs): + """Build a registered incoherent model around a coherent crystal. + + Direct construction remains the ordinary Python API; this exists so a + stored type name can be turned into a class without importing an + arbitrary path from a file. + + :param str key: + Registry key. + :param CTRcalc.SXRDCrystal crystal: + Coherent crystal to wrap. + :param keyargs: + Forwarded to the model constructor. + :rtype: IncoherentModel + """ + return _registered(key).model_class(crystal, **keyargs) + + +def create_incoherent_model_from_config(crystal, config): + """Rebuild a model from :meth:`IncoherentModel.to_config`. + + Only registered keys can be instantiated, and only wrapper-local state is + restored: the coherent crystal is supplied by the caller and keeps the + parameters it already has. + + :param CTRcalc.SXRDCrystal crystal: + Coherent crystal to wrap. + :param dict config: + Dictionary produced by ``to_config``. + :rtype: IncoherentModel + :raises ValueError: + If the type is unknown or the dictionary carries unexpected keys. + """ + if not isinstance(config, dict): + raise ValueError("An incoherent model configuration must be a dict") + unexpected = set(config) - {"type", "settings", "parameters"} + if unexpected: + raise ValueError( + "Unexpected keys in incoherent model configuration: " + + ", ".join(sorted(unexpected)) + ) + if "type" not in config: + raise ValueError("An incoherent model configuration needs a type") + + record = _registered(config["type"]) + settings = config.get("settings") or {} + model = record.model_class(crystal, **dict(settings)) + parameters = config.get("parameters") + if parameters: + model._local_parameters_from_dict(parameters) + return model diff --git a/orgui/datautils/xrayutils/test/test_CTRincoherent.py b/orgui/datautils/xrayutils/test/test_CTRincoherent.py new file mode 100644 index 0000000..68cb067 --- /dev/null +++ b/orgui/datautils/xrayutils/test/test_CTRincoherent.py @@ -0,0 +1,667 @@ +"""Contract tests for the incoherent CTR wrapper. + +Every model here is synthetic. Increment 3 of +``doc/design/incoherent_ctr_models.md`` requires the contract to hold without +the Poisson implementation, so nothing in this file constructs or imports a +``PoissonSurface``. +""" + +import copy +import unittest + +import numpy as np + +from .. import CTRcalc, CTRincoherent +from ..CTRincoherent import ( + CoherentStateEnsembleModel, + IncoherentF2Model, + IncoherentState, + available_incoherent_models, + create_incoherent_model, + create_incoherent_model_from_config, + register_incoherent_model, + unregister_incoherent_model, +) +from . import _poisson_oracle + + +def _crystal(fitted=0): + """Return a plain coherent crystal with ``fitted`` crystal parameters.""" + film = CTRcalc.UnitCell([3.0, 3.0, 6.0], [90.0, 90.0, 90.0], name="film") + film.addAtom("C", [0.0, 0.0, 0.0], 0.1, 0.1, 1.0, layer=0) + crystal = CTRcalc.SXRDCrystal( + _poisson_oracle.layered_cell(1, "bulk"), film, stacking=np.array([1]) + ) + # Crystal fit parameters come from the component cell, which is how a + # real crystal acquires its tail. + definitions = (("iDW", (0.0, 5.0)), ("occ", (0.0, 1.0))) + for index in range(fitted): + parameter, limits = definitions[index] + crystal["film"].addFitParameter( + (0, parameter), limits=limits, name=f"film {parameter}" + ) + return crystal + + +class _NoLocalModel(IncoherentF2Model): + """Wrapper with no local fit parameters.""" + + model_type = "test_no_local" + + def _evaluate_F2(self, context): # noqa: N802 + return np.abs(context.coherent.total) ** 2 + + +class _OneLocalModel(IncoherentF2Model): + """Wrapper with one local fit parameter.""" + + model_type = "test_one_local" + parameterLookup = {"scale": 0} + parameterLookup_inv = {0: "scale"} + + def __init__(self, crystal, *, scale=0.5, name="one"): + super().__init__(crystal, name=name) + self.basis = np.array([float(scale)]) + self.basis_0 = np.array([float(scale)]) + self.errors = None + + def _config_settings(self): + return {} + + def _evaluate_F2(self, context): # noqa: N802 + return self.basis[0] * np.abs(context.coherent.total) ** 2 + + +class _ThreeLocalModel(IncoherentF2Model): + """Wrapper with three local fit parameters.""" + + model_type = "test_three_local" + parameterLookup = {"a": 0, "b": 1, "c": 2} + parameterLookup_inv = {0: "a", 1: "b", 2: "c"} + + def __init__(self, crystal, *, name="three"): + super().__init__(crystal, name=name) + self.basis = np.array([0.1, 0.2, 0.3]) + self.basis_0 = np.copy(self.basis) + self.errors = None + + def _evaluate_F2(self, context): # noqa: N802 + return float(np.sum(self.basis)) * np.abs(context.coherent.total) ** 2 + + +class _TwoStateSource: + """State source standing in for a real surface decomposition.""" + + def __init__(self, amplitudes, probabilities): + self._amplitudes = amplitudes + self._probabilities = probabilities + + def iter_states(self): + """Yield the states one at a time.""" + for index, (amplitude, probability) in enumerate( + zip(self._amplitudes, self._probabilities) + ): + yield IncoherentState(index, probability, amplitude) + + +class _TwoStateEnsemble(CoherentStateEnsembleModel): + """Ensemble of two explicitly supplied complete state amplitudes.""" + + model_type = "test_two_state" + + def __init__(self, crystal, *, kappa=0.0, offsets=(5.0, -5.0), name="two"): + super().__init__(crystal, name=name) + self.basis = np.array([float(kappa)]) + self.basis_0 = np.array([float(kappa)]) + self.errors = None + self.offsets = tuple(offsets) + self.probabilities = (0.5, 0.5) + + parameterLookup = {"incoherent_fraction": 0} + parameterLookup_inv = {0: "incoherent_fraction"} + + @property + def incoherent_fraction(self): + """Return the mixing fraction held in the local basis.""" + return float(self.basis[0]) + + def state_amplitudes(self, context): + """Return the two complete state amplitudes for this context.""" + total = context.coherent.total + return [total + offset for offset in self.offsets] + + def _iter_states(self, context): + source = _TwoStateSource( + self.state_amplitudes(context), self.probabilities + ) + return source.iter_states() + + +class ContractMixin(unittest.TestCase): + """Shared coordinates for the contract tests.""" + + H = np.zeros(4) + K = np.zeros(4) + L = np.array([0.4, 0.9, 1.5, 2.2]) + + +class TestParameterBlockContract(ContractMixin): + """The composite fit API keeps one stable local-then-coherent order.""" + + def models(self): + """Yield one wrapper of each local-parameter count.""" + no_local = _NoLocalModel(_crystal(fitted=2)) + + one_local = _OneLocalModel(_crystal(fitted=2)) + one_local.addFitParameter("scale", limits=(0.0, 1.0)) + + three_local = _ThreeLocalModel(_crystal(fitted=1)) + for parameter in ("a", "b", "c"): + three_local.addFitParameter(parameter, limits=(0.0, 1.0)) + + return { + "zero local": (no_local, 0, 2), + "one local": (one_local, 1, 2), + "three local": (three_local, 3, 1), + } + + def test_every_contract_array_has_the_same_length_and_order(self): + """One check catches any composite override left un-overridden. + + ``LinearFitFunctions`` implements these as local-only methods, so a + forgotten override does not raise: it returns a vector one block + short. Comparing every array at once is what makes that visible. + """ + for label, (model, n_local, n_coherent) in self.models().items(): + with self.subTest(case=label): + total = n_local + n_coherent + start, lower, upper = model.getStartParamAndLimits() + self.assertEqual(len(model.fitparnames), total) + self.assertEqual(len(model.getInitialParameters()), total) + self.assertEqual(len(model.parameter_list()), total) + self.assertEqual(len(start), total) + self.assertEqual(len(lower), total) + self.assertEqual(len(upper), total) + self.assertEqual(model.n_local_parameters, n_local) + self.assertEqual(model.n_coherent_parameters, n_coherent) + + # The coherent block is the tail, in the crystal's own order. + crystal_names = list(model.coherent_model.fitparnames) + self.assertEqual( + model.fitparnames[n_local:], crystal_names + ) + np.testing.assert_allclose( + model.getInitialParameters()[n_local:], + model.coherent_model.getInitialParameters(), + ) + + def test_priors_cover_both_blocks(self): + """Priors follow the same order and length as the names.""" + model = _OneLocalModel(_crystal(fitted=2)) + model.addFitParameter("scale", limits=(0.0, 1.0)) + self.assertEqual(len(model.priors), len(model.fitparnames)) + # Prior entries may be arrays, so compare them element-wise. + for mine, theirs in zip( + model.priors[1:], list(model.coherent_model.priors) + ): + np.testing.assert_array_equal(mine, theirs) + np.testing.assert_array_equal(model.priors[0], (0.0, 1.0)) + + def test_set_parameters_splits_the_vector(self): + """Each block receives its own segment.""" + model = _OneLocalModel(_crystal(fitted=2)) + model.addFitParameter("scale", limits=(0.0, 1.0)) + + model.setParameters([0.75, 1.5, 2.5]) + self.assertAlmostEqual(model.basis[0], 0.75) + np.testing.assert_allclose( + model.coherent_model.getInitialParameters(), [1.5, 2.5] + ) + np.testing.assert_allclose( + model.getInitialParameters(), [0.75, 1.5, 2.5] + ) + + def test_bad_lengths_leave_the_model_untouched(self): + """The whole vector is validated before either block is written.""" + model = _OneLocalModel(_crystal(fitted=2)) + model.addFitParameter("scale", limits=(0.0, 1.0)) + before = model.getInitialParameters().copy() + + for bad in ([0.5], [0.5, 1.0], [0.5, 1.0, 2.0, 3.0]): + with self.subTest(length=len(bad)): + with self.assertRaisesRegex(ValueError, "entries but this model"): + model.setParameters(bad) + np.testing.assert_allclose( + model.getInitialParameters(), before + ) + + def test_local_only_setter_refuses(self): + """``setFitParameters`` would silently drop the coherent tail.""" + model = _OneLocalModel(_crystal(fitted=2)) + model.addFitParameter("scale", limits=(0.0, 1.0)) + with self.assertRaisesRegex(NotImplementedError, "setParameters"): + model.setFitParameters([0.5]) + + def test_errors_round_trip_across_both_blocks(self): + """Errors split like parameters, and clear together.""" + model = _OneLocalModel(_crystal(fitted=2)) + model.addFitParameter("scale", limits=(0.0, 1.0)) + + model.setFitErrors([0.01, 0.02, 0.03]) + np.testing.assert_allclose( + model.getFitErrors(), [0.01, 0.02, 0.03] + ) + + model.setFitErrors(None) + with self.assertRaises(ValueError): + model.getFitErrors() + + def test_empty_local_block_is_skipped_not_queried(self): + """A block with no parameters must not raise for missing errors.""" + model = _NoLocalModel(_crystal(fitted=2)) + model.setFitErrors([0.05, 0.06]) + np.testing.assert_allclose(model.getFitErrors(), [0.05, 0.06]) + + def test_bounds_reach_both_blocks(self): + """``setLimits`` splits like the parameter vector.""" + model = _OneLocalModel(_crystal(fitted=2)) + model.addFitParameter("scale", limits=(0.0, 1.0)) + model.setLimits([(0.1, 0.9), (0.2, 4.0), (0.3, 4.5)]) + _, lower, upper = model.getStartParamAndLimits() + np.testing.assert_allclose(lower, [0.1, 0.2, 0.3]) + np.testing.assert_allclose(upper, [0.9, 4.0, 4.5]) + + def test_wrapping_requires_a_crystal(self): + """The wrapper owns one primary coherent model.""" + with self.assertRaisesRegex(TypeError, "wraps one SXRDCrystal"): + _NoLocalModel(object()) + + +class TestSerializationAndRegistry(ContractMixin): + """Copying, configuration round-trips, and registry behavior.""" + + def test_deepcopy_keeps_both_blocks(self): + """``CTROptimizer`` deep-copies the model once per fit.""" + model = _OneLocalModel(_crystal(fitted=2)) + model.addFitParameter("scale", limits=(0.0, 1.0)) + model.setParameters([0.75, 1.5, 2.5]) + + copied = copy.deepcopy(model) + np.testing.assert_allclose( + copied.getInitialParameters(), model.getInitialParameters() + ) + self.assertEqual(copied.fitparnames, model.fitparnames) + + copied.setParameters([0.25, 1.0, 2.0]) + np.testing.assert_allclose( + model.getInitialParameters(), [0.75, 1.5, 2.5] + ) + + def test_parameter_dict_round_trip(self): + """Both subtrees survive a dictionary round-trip.""" + model = _OneLocalModel(_crystal(fitted=2)) + model.addFitParameter("scale", limits=(0.0, 1.0)) + model.setParameters([0.75, 1.5, 2.5]) + stored = model.parametersToDict() + self.assertEqual(set(stored), {"local", "coherent"}) + + restored = _OneLocalModel(_crystal(fitted=2)) + restored.parametersFromDict(stored) + self.assertEqual(restored.fitparnames, model.fitparnames) + np.testing.assert_allclose( + restored.getInitialParameters(), model.getInitialParameters() + ) + + def test_clear_parameters_clears_both(self): + """Clearing is composite, like the other parameter operations.""" + model = _OneLocalModel(_crystal(fitted=2)) + model.addFitParameter("scale", limits=(0.0, 1.0)) + model.clearParameters() + self.assertEqual(model.fitparnames, []) + self.assertEqual(model.n_local_parameters, 0) + self.assertEqual(model.n_coherent_parameters, 0) + + def test_config_round_trip_restores_only_local_state(self): + """The factory wraps a supplied crystal and restores local state.""" + register_incoherent_model(_OneLocalModel) + try: + model = _OneLocalModel(_crystal(fitted=2)) + model.addFitParameter("scale", limits=(0.0, 1.0)) + model.setParameters([0.75, 1.5, 2.5]) + config = model.to_config() + self.assertEqual(config["type"], "test_one_local") + self.assertEqual(set(config), {"type", "settings", "parameters"}) + + crystal = _crystal(fitted=2) + restored = create_incoherent_model_from_config(crystal, config) + self.assertIs(restored.coherent_model, crystal) + self.assertEqual( + restored.fitparnames[: restored.n_local_parameters], + ["one scale"], + ) + self.assertAlmostEqual(restored.basis[0], 0.75) + finally: + unregister_incoherent_model("test_one_local") + + def test_registry_rejects_bad_entries(self): + """Only concrete, keyed, quantity-bearing subclasses register.""" + register_incoherent_model(_NoLocalModel) + try: + self.assertIn("test_no_local", available_incoherent_models()) + info = available_incoherent_models()["test_no_local"] + self.assertIs(info.model_class, _NoLocalModel) + self.assertEqual(info.output_quantity, "F2") + self.assertIn("kinematical", info.supported_forward_models) + + with self.assertRaisesRegex(ValueError, "already registered"): + register_incoherent_model(_NoLocalModel) + + with self.assertRaisesRegex(TypeError, "IncoherentModel"): + register_incoherent_model(dict) + + with self.assertRaisesRegex(TypeError, "abstract"): + register_incoherent_model(IncoherentF2Model) + + class _Unkeyed(_NoLocalModel): + model_type = "" + + with self.assertRaisesRegex(ValueError, "model_type"): + register_incoherent_model(_Unkeyed) + finally: + unregister_incoherent_model("test_no_local") + + def test_registry_mapping_is_read_only(self): + """Discovery returns a view, not the live registry.""" + models = available_incoherent_models() + with self.assertRaises(TypeError): + models["anything"] = None + + def test_unknown_keys_fail_deterministically(self): + """Deserialization instantiates registered keys only.""" + with self.assertRaisesRegex(ValueError, "Unknown incoherent model"): + create_incoherent_model("not_registered", _crystal()) + with self.assertRaisesRegex(ValueError, "Unknown incoherent model"): + create_incoherent_model_from_config( + _crystal(), {"type": "not_registered"} + ) + with self.assertRaisesRegex(ValueError, "Unexpected keys"): + create_incoherent_model_from_config( + _crystal(), {"type": "x", "surprise": 1} + ) + + def test_third_party_models_need_no_core_changes(self): + """A model defined outside the package registers and evaluates.""" + + @register_incoherent_model + class _ThirdParty(IncoherentF2Model): + """Externally defined model.""" + + model_type = "test_third_party" + + def _evaluate_F2(self, context): # noqa: N802 + return 2.0 * np.abs(context.coherent.total) ** 2 + + try: + model = create_incoherent_model("test_third_party", _crystal()) + np.testing.assert_allclose( + model.F2(self.H, self.K, self.L), + 2.0 * model.coherent_model.F2(self.H, self.K, self.L), + rtol=1e-14, + ) + finally: + unregister_incoherent_model("test_third_party") + + +class TestF2Boundary(ContractMixin): + """The squared-structure-factor boundary and its guarantees.""" + + def test_F2_is_real_and_nonnegative(self): + """The public boundary validates what it returns.""" + model = _NoLocalModel(_crystal()) + value = model.F2(self.H, self.K, self.L) + self.assertEqual(value.dtype, np.float64) + self.assertTrue(np.all(value >= 0.0)) + np.testing.assert_allclose( + value, model.coherent_model.F2(self.H, self.K, self.L), rtol=1e-14 + ) + + def test_no_F_method_is_exposed(self): + """A mixed state has no unique complex amplitude.""" + model = _NoLocalModel(_crystal()) + self.assertFalse(hasattr(model, "F")) + + def test_negative_or_nonfinite_results_are_rejected(self): + """A broken model fails at the boundary rather than downstream.""" + + class _Negative(IncoherentF2Model): + """Deliberately invalid model.""" + + model_type = "test_negative" + + def _evaluate_F2(self, context): # noqa: N802 + return -np.ones_like(np.asarray(context.coordinates[2])) + + class _NotFinite(IncoherentF2Model): + """Deliberately invalid model.""" + + model_type = "test_not_finite" + + def _evaluate_F2(self, context): # noqa: N802 + return np.full_like( + np.asarray(context.coordinates[2]), np.nan + ) + + with self.assertRaisesRegex(ValueError, "negative F2"): + _Negative(_crystal()).F2(self.H, self.K, self.L) + with self.assertRaisesRegex(ValueError, "non-finite F2"): + _NotFinite(_crystal()).F2(self.H, self.K, self.L) + + def test_validate_rejects_unsupported_forward_models(self): + """An F2 wrapper must not be reinterpreted as a reflectivity.""" + model = _NoLocalModel(_crystal()) + model.validate(forward_model="kinematical") + with self.assertRaisesRegex(ValueError, "does not support"): + model.validate(forward_model="dwba") + + +class TestEnsembleMixing(ContractMixin): + """Endpoints and interpolation of the generic state ensemble.""" + + def ensemble(self, kappa, offsets=(5.0, -5.0)): + """Return a two-state ensemble with the given mixing fraction.""" + return _TwoStateEnsemble(_crystal(), kappa=kappa, offsets=offsets) + + def test_coherent_endpoint_matches_the_wrapped_crystal(self): + """``kappa = 0`` is the live coherent result, exactly.""" + model = self.ensemble(0.0) + np.testing.assert_allclose( + model.F2(self.H, self.K, self.L), + model.coherent_model.F2(self.H, self.K, self.L), + rtol=1e-14, + ) + + def test_incoherent_endpoint_is_the_squared_average(self): + """``kappa = 1`` averages the squared state amplitudes.""" + model = self.ensemble(1.0) + context = CTRincoherent.KinematicIncoherentContext( + model.coherent_model, self.H, self.K, self.L + ) + amplitudes = model.state_amplitudes(context) + expected = sum( + 0.5 * np.abs(amplitude) ** 2 for amplitude in amplitudes + ) + np.testing.assert_allclose( + model.F2(self.H, self.K, self.L), expected, rtol=1e-12 + ) + + def test_partial_mixing_is_the_convex_interpolation(self): + """Intermediate fractions interpolate the two endpoints linearly.""" + coherent = self.ensemble(0.0).F2(self.H, self.K, self.L) + incoherent = self.ensemble(1.0).F2(self.H, self.K, self.L) + for kappa in (0.25, 0.5, 0.75): + with self.subTest(kappa=kappa): + np.testing.assert_allclose( + self.ensemble(kappa).F2(self.H, self.K, self.L), + (1.0 - kappa) * coherent + kappa * incoherent, + rtol=1e-12, + ) + + def test_the_endpoints_actually_differ(self): + """Otherwise the interpolation test would pass on a constant.""" + coherent = self.ensemble(0.0).F2(self.H, self.K, self.L) + incoherent = self.ensemble(1.0).F2(self.H, self.K, self.L) + self.assertGreater( + np.max(np.abs(incoherent - coherent)) / np.max(incoherent), 0.01 + ) + # The incoherent endpoint drops a cross term, so it never falls below. + self.assertTrue(np.all(incoherent >= coherent - 1e-9)) + + def test_a_single_state_is_independent_of_kappa(self): + """With one populated height there is no variance to expose.""" + + class _SingleState(_TwoStateEnsemble): + """One-state ensemble.""" + + model_type = "test_single_state" + + def __init__(self, crystal, *, kappa=0.0): + super().__init__(crystal, kappa=kappa, offsets=(0.0,)) + self.probabilities = (1.0,) + + values = [ + _SingleState(_crystal(), kappa=kappa).F2(self.H, self.K, self.L) + for kappa in (0.0, 0.5, 1.0) + ] + for value in values[1:]: + np.testing.assert_allclose(value, values[0], rtol=1e-12) + + def test_unnormalized_masses_are_rejected(self): + """Probabilities must be normalized over the retained states once.""" + + class _Unnormalized(_TwoStateEnsemble): + """Ensemble whose masses do not sum to one.""" + + model_type = "test_unnormalized" + + def __init__(self, crystal): + super().__init__(crystal, kappa=1.0) + self.probabilities = (0.5, 0.2) + + with self.assertRaisesRegex(ValueError, "sum to"): + _Unnormalized(_crystal()).F2(self.H, self.K, self.L) + + def test_invalid_mixing_fractions_are_rejected(self): + """``kappa`` is dimensionless and lives in ``[0, 1]``.""" + for kappa in (-0.1, 1.5, np.nan): + with self.subTest(kappa=kappa): + model = self.ensemble(1.0) + model.basis[0] = kappa + with self.assertRaisesRegex(ValueError, "incoherent_fraction"): + model.F2(self.H, self.K, self.L) + + def test_reconstruction_diagnostic_is_available(self): + """The state sum is a diagnostic against the live coherent result.""" + model = self.ensemble(1.0, offsets=(2.0, -2.0)) + self.assertAlmostEqual( + model.coherence_reconstruction_error(self.H, self.K, self.L), + 0.0, + places=12, + ) + + +class TestEvaluationContext(ContractMixin): + """The lazy decomposition shared by every model.""" + + def test_the_coherent_decomposition_is_evaluated_once(self): + """Streaming many states must not re-evaluate the common crystal.""" + crystal = _crystal() + calls = [] + original = crystal.evaluate_kinematic + + def counting(*args, **keyargs): + calls.append(args) + return original(*args, **keyargs) + + crystal.evaluate_kinematic = counting + context = CTRincoherent.KinematicIncoherentContext( + crystal, self.H, self.K, self.L + ) + for _ in range(5): + context.coherent # noqa: B018 + context.component("film") + context.common_amplitude("film") + self.assertEqual(len(calls), 1) + + def test_common_amplitude_removes_only_the_named_component(self): + """The common part is everything except the target component.""" + crystal = _crystal() + context = CTRincoherent.KinematicIncoherentContext( + crystal, self.H, self.K, self.L + ) + np.testing.assert_allclose( + context.common_amplitude("film"), context.coherent.bulk, rtol=1e-14 + ) + + def test_states_are_evaluated_at_transformed_coordinates(self): + """An outer domain changes where a component state must be evaluated. + + The context must call the evaluator once per domain at that domain's + coordinates, and combine the results with the crystal's own area, + weight, and occupancy scaling. + """ + crystal = _crystal() + transforms = [ + (np.identity(3), 0.6), + (np.diag([1.0, 1.0, 0.5]), 0.4), + ] + crystal.setDomain(0, transforms) + context = CTRincoherent.KinematicIncoherentContext( + crystal, self.H, self.K, self.L + ) + + seen = [] + + def evaluator(component, h, k, l): # noqa: E741 + seen.append(np.copy(np.asarray(l))) + return _TwoStateSource( + [ + np.asarray(l, dtype=np.complex128), + np.zeros_like(l, dtype=np.complex128), + ], + [0.5, 0.5], + ) + + states = list(context.iter_component_states("film", evaluator)) + + self.assertEqual(len(seen), 2) + np.testing.assert_allclose(seen[0], self.L) + np.testing.assert_allclose(seen[1], 0.5 * self.L) + + common = context.common_amplitude("film") + area_weight = crystal.reference_area / crystal.uc_surface_list[0].uc_area + expected = common + area_weight * ( + 0.6 * self.L + 0.4 * (0.5 * self.L) + ) + self.assertEqual(len(states), 2) + np.testing.assert_allclose(states[0].amplitude, expected, rtol=1e-12) + self.assertEqual(states[0].probability, 0.5) + + def test_inconsistent_state_metadata_is_rejected(self): + """Only the amplitudes may depend on the domain transform.""" + crystal = _crystal() + crystal.setDomain(0, [(np.identity(3), 0.5), (np.identity(3), 0.5)]) + context = CTRincoherent.KinematicIncoherentContext( + crystal, self.H, self.K, self.L + ) + masses = iter([(0.5, 0.5), (0.25, 0.75)]) + + def evaluator(component, h, k, l): # noqa: E741 + zeros = np.zeros_like(np.asarray(l), dtype=np.complex128) + return _TwoStateSource([zeros, zeros], next(masses)) + + with self.assertRaisesRegex(ValueError, "disagree on state"): + list(context.iter_component_states("film", evaluator)) + + +if __name__ == "__main__": + unittest.main() From 0391b2f0166975581f99f4c001b63771aaf0bb5e Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Fri, 11 Sep 2026 19:46:11 -0400 Subject: [PATCH 10/13] feat(phys): implement the Poisson height-domain model Mixes the complete crystal amplitude of each retained flat height, so bulk-surface interference stays inside every domain. The coherent endpoint is the cached coherent result, not a state-sum reconstruction. --- orgui/datautils/xrayutils/CTRincoherent.py | 172 ++++++++ .../xrayutils/test/test_CTRincoherent.py | 366 +++++++++++++++++- 2 files changed, 534 insertions(+), 4 deletions(-) diff --git a/orgui/datautils/xrayutils/CTRincoherent.py b/orgui/datautils/xrayutils/CTRincoherent.py index c393368..9e5a332 100644 --- a/orgui/datautils/xrayutils/CTRincoherent.py +++ b/orgui/datautils/xrayutils/CTRincoherent.py @@ -55,6 +55,7 @@ import numpy as np from .CTRcalc import SXRDCrystal +from .CTRfilm import PoissonSurface from .CTRutil import LinearFitFunctions __all__ = [ @@ -64,6 +65,7 @@ "IncoherentModelInfo", "IncoherentState", "KinematicIncoherentContext", + "PoissonHeightDomains", "available_incoherent_models", "create_incoherent_model", "create_incoherent_model_from_config", @@ -834,3 +836,173 @@ def create_incoherent_model_from_config(crystal, config): if parameters: model._local_parameters_from_dict(parameters) return model + + +@register_incoherent_model(builtin=True) +class PoissonHeightDomains(CoherentStateEnsembleModel): + """Laterally large height domains on a Poisson-roughened surface. + + Each coherence patch sees one flat height drawn from the surface's + ``PoissonProfile``. The states are the complete coherent amplitudes of the + crystal at those heights, so bulk-surface and Film-surface interference + stays inside each domain, and ``incoherent_fraction`` interpolates between + the coherent limit and their squared average. + + Constructing this wrapper is the opt-in: the ``PoissonSurface`` it targets + stays an ordinary coherent structural model on its own. + + :param CTRcalc.SXRDCrystal crystal: + Coherent crystal to wrap. + :param str surface: + Stable name of the target ``PoissonSurface`` component. + :param float incoherent_fraction: + Dimensionless mixing fraction ``kappa`` in ``[0, 1]``. Zero is fully + coherent, one is the large-domain limit. + :param int exact_layer_count: + Retain every height state when the profile populates no more than + this many. A policy switch, not a cap. + :param str name: + Wrapper name, the default prefix for its local parameter. + :raises ValueError: + If the target is missing, duplicated, not a ``PoissonSurface``, not + the topmost component, or the fraction is out of range. + """ + + model_type = "poisson_height_domains" + parameterLookup = {"incoherent_fraction": 0} + parameterLookup_inv = {0: "incoherent_fraction"} + + def __init__( + self, + crystal, + *, + surface, + incoherent_fraction=0.0, + exact_layer_count=10, + name="incoherent", + ): + super().__init__(crystal, name=name) + self.surface = str(surface) + if ( + int(exact_layer_count) != exact_layer_count + or exact_layer_count < 1 + ): + raise ValueError("exact_layer_count must be a positive integer") + self.exact_layer_count = int(exact_layer_count) + # One local basis entry, so a fixed fraction and a fitted one are the + # same stored value. `basis_0` is seeded too, because + # `updateFromParameters` rebuilds `basis` from it whenever the + # fraction is not itself a fit parameter. + self.basis = np.array([0.0]) + self.basis_0 = np.array([0.0]) + self.errors = None + self.incoherent_fraction = incoherent_fraction + self._validate_target() + + @property + def incoherent_fraction(self): + """Return the dimensionless mixing fraction held in the basis.""" + return float(self.basis[0]) + + @incoherent_fraction.setter + def incoherent_fraction(self, value): + """Set the mixing fraction, which must be finite and in ``[0, 1]``.""" + value = float(value) + if not np.isfinite(value) or not 0.0 <= value <= 1.0: + raise ValueError( + f"incoherent_fraction must be finite and within [0, 1], " + f"got {value}" + ) + self.basis[0] = value + self.basis_0[0] = value + + def target_surface(self): + """Return the targeted ``PoissonSurface`` component. + + Selection is by stable name: component identities do not survive the + deep copy the optimizer performs once per fit. + + :rtype: CTRfilm.PoissonSurface + :raises ValueError: + If no component, or more than one, carries the name. + """ + components = self._coherent_model.uc_surface_list + matches = [ + component + for component in components + if getattr(component, "name", None) == self.surface + ] + if not matches: + available = ", ".join( + str(getattr(component, "name", "?")) + for component in components + ) + raise ValueError( + f"No component named {self.surface!r} on the crystal. " + f"Available: {available}" + ) + if len(matches) > 1: + raise ValueError( + f"{len(matches)} components are named {self.surface!r}; " + "the target must be identified by a unique name" + ) + return matches[0] + + def _validate_target(self): + """Check the target type and the supported stacking topology.""" + target = self.target_surface() + if not isinstance(target, PoissonSurface): + raise ValueError( + f"Component {self.surface!r} is a " + f"{type(target).__name__}, not a PoissonSurface" + ) + + crystal = self._coherent_model + ordered = list(crystal.uc_surface_list_ordered) + if not ordered or ordered[-1] is not target: + raise ValueError( + f"{self.surface!r} must be the topmost component in stacking " + "order. A component above it could have a position which " + "depends on the selected height, so it cannot be treated as " + "part of the common amplitude" + ) + levels = np.asarray(crystal.uc_stacking_ordered) + if levels.size and int(np.sum(levels == levels[-1])) != 1: + raise ValueError( + f"{self.surface!r} shares its stacking level with another " + "component, whose position could depend on the selected " + "height" + ) + + def _validate_model(self, forward_model): + """Check that the target is bound to the Film it corrects.""" + target = self._bound_target() + if target.underlying_film is None: + raise ValueError( + f"{self.surface!r} is not stacked immediately above a Film, " + "so it has no sharp boundary to correct" + ) + + def _bound_target(self): + """Return the target after applying any pending stacking.""" + crystal = self._coherent_model + if crystal.enable_uc_stacking: + crystal.apply_stacking() + return self.target_surface() + + def _iter_states(self, context): + """Stream complete crystal amplitudes, one flat height at a time.""" + + def state_evaluator(component, h, k, l): # noqa: E741 + return component.flat_domain_corrections( + h, k, l, exact_layer_count=self.exact_layer_count + ) + + return context.iter_component_states(self.surface, state_evaluator) + + def _config_settings(self): + """Return the non-basis settings. The fraction is not among them.""" + return { + "surface": self.surface, + "exact_layer_count": self.exact_layer_count, + } diff --git a/orgui/datautils/xrayutils/test/test_CTRincoherent.py b/orgui/datautils/xrayutils/test/test_CTRincoherent.py index 68cb067..ad1f9e3 100644 --- a/orgui/datautils/xrayutils/test/test_CTRincoherent.py +++ b/orgui/datautils/xrayutils/test/test_CTRincoherent.py @@ -1,9 +1,11 @@ -"""Contract tests for the incoherent CTR wrapper. +"""Tests for the incoherent CTR wrapper and the Poisson height-domain model. -Every model here is synthetic. Increment 3 of +The contract tests use synthetic models only. Increment 3 of ``doc/design/incoherent_ctr_models.md`` requires the contract to hold without -the Poisson implementation, so nothing in this file constructs or imports a -``PoissonSurface``. +the Poisson implementation, so nothing above ``TestPoissonHeightDomains`` +constructs a ``PoissonSurface``. The increment-4 classes at the end of the +file do, and check the model against the independently stacked flat-height +crystals of ``_poisson_oracle``. """ import copy @@ -22,6 +24,7 @@ register_incoherent_model, unregister_incoherent_model, ) +from ..CTRdistributions import PoissonProfile from . import _poisson_oracle @@ -663,5 +666,360 @@ def evaluator(component, h, k, l): # noqa: E741 list(context.iter_component_states("film", evaluator)) +class PoissonModelMixin(ContractMixin): + """Builders for the Poisson height-domain model.""" + + def model(self, profile, kappa=0.0, n_layers=2, **keyargs): + """Return the model, its crystal, its surface, and the base width.""" + w_base = _poisson_oracle.minimum_base_width(profile) + crystal, surface = _poisson_oracle.poisson_crystal( + profile, w_base=w_base, n_layers=n_layers, **keyargs + ) + crystal.apply_stacking() + model = CTRincoherent.PoissonHeightDomains( + crystal, surface="surface", incoherent_fraction=kappa + ) + return model, crystal, surface, w_base + + def oracle_incoherent(self, surface, w_base, n_layers=2, **keyargs): + """Return the state average of independently stacked crystals.""" + states = surface.flat_domain_corrections(self.H, self.K, self.L) + total = None + for layer, probability in zip( + states.layer_numbers, states.probabilities + ): + amplitude = _poisson_oracle.flat_height_crystal( + int(layer), w_base=w_base, n_layers=n_layers, **keyargs + ).F(self.H, self.K, self.L) + squared = probability * np.abs(amplitude) ** 2 + total = squared if total is None else total + squared + return total + + +class TestPoissonHeightDomains(PoissonModelMixin): + """Endpoints and interpolation against independent flat-height crystals.""" + + PROFILES = { + "poisson growth": PoissonProfile( + 1.5, alpha=0.6, tail_probability=1e-14 + ), + "poisson etching": PoissonProfile( + -1.5, alpha=0.6, tail_probability=1e-14 + ), + "fractional step": PoissonProfile( + 1.5, alpha=0.0, tail_probability=1e-14 + ), + "step and poisson": PoissonProfile( + 2.3, alpha=0.4, offset=0.4, tail_probability=1e-14 + ), + } + + def test_coherent_endpoint_is_the_wrapped_crystal(self): + """``kappa = 0`` must not be perturbed by the finite state support. + + It comes from the coherent result already cached in the context, so + it is the crystal's own ``F2`` rather than a state-sum + reconstruction. + """ + for label, profile in self.PROFILES.items(): + with self.subTest(case=label): + model, crystal, _, _ = self.model(profile, kappa=0.0) + np.testing.assert_array_equal( + model.F2(self.H, self.K, self.L), + crystal.F2(self.H, self.K, self.L), + ) + + def test_incoherent_endpoint_matches_deterministic_crystals(self): + """``kappa = 1`` is the weighted ``F2`` average of flat crystals.""" + for label, profile in self.PROFILES.items(): + with self.subTest(case=label): + model, _, surface, w_base = self.model(profile, kappa=1.0) + np.testing.assert_allclose( + model.F2(self.H, self.K, self.L), + self.oracle_incoherent(surface, w_base), + rtol=1e-12, + ) + + def test_partial_mixing_is_the_convex_interpolation(self): + """Intermediate fractions fill in linearly between the endpoints.""" + profile = self.PROFILES["poisson growth"] + coherent = self.model(profile, kappa=0.0)[0].F2(self.H, self.K, self.L) + incoherent = self.model(profile, kappa=1.0)[0].F2( + self.H, self.K, self.L + ) + self.assertGreater( + np.max(np.abs(incoherent - coherent)) / np.max(incoherent), 0.05 + ) + for kappa in (0.2, 0.35, 0.8): + with self.subTest(kappa=kappa): + model, _, _, _ = self.model(profile, kappa=kappa) + np.testing.assert_allclose( + model.F2(self.H, self.K, self.L), + (1.0 - kappa) * coherent + kappa * incoherent, + rtol=1e-12, + ) + + def test_the_anti_bragg_minimum_fills_in_with_kappa(self): + """Two equally populated heights cancel coherently, not squared. + + ``PoissonProfile(0.5, alpha=0.0)`` populates layers -1 and 0 at + exactly one half each, so the coherent minimum along the rod is the + deepest the interpolation has to fill. + """ + profile = PoissonProfile(0.5, alpha=0.0, tail_probability=1e-14) + scan = np.linspace(0.05, 3.0, 200) + zeros = np.zeros_like(scan) + + w_base = 4.0 + crystal, surface = _poisson_oracle.poisson_crystal( + profile, w_base=w_base + ) + crystal.apply_stacking() + values = {} + for kappa in (0.0, 0.5, 1.0): + model = CTRincoherent.PoissonHeightDomains( + crystal, surface="surface", incoherent_fraction=kappa + ) + values[kappa] = model.F2(zeros, zeros, scan) + + # The minimum is filled in, and half way there it is exactly half way. + deepest = int(np.argmin(values[0.0])) + self.assertGreater(values[1.0][deepest], values[0.0][deepest]) + np.testing.assert_allclose( + values[0.5], 0.5 * values[0.0] + 0.5 * values[1.0], rtol=1e-12 + ) + + def test_a_deterministic_height_is_independent_of_kappa(self): + """One populated state leaves no variance for the mixture to expose.""" + profile = PoissonProfile(2.0, alpha=0.0, tail_probability=1e-14) + reference = None + for kappa in (0.0, 0.4, 1.0): + with self.subTest(kappa=kappa): + model, _, surface, _ = self.model(profile, kappa=kappa) + states = surface.flat_domain_corrections( + self.H, self.K, self.L + ) + self.assertEqual(states.layer_numbers.size, 1) + value = model.F2(self.H, self.K, self.L) + if reference is None: + reference = value + np.testing.assert_allclose(value, reference, rtol=1e-12) + + def test_bulk_interference_stays_inside_each_state(self): + """Averaging only the surface correction is a different quantity.""" + profile = PoissonProfile(0.5, alpha=0.0, tail_probability=1e-14) + model, _, surface, w_base = self.model(profile, kappa=1.0) + states = surface.flat_domain_corrections(self.H, self.K, self.L) + common = _poisson_oracle.flat_height_crystal(-1, w_base=w_base).F( + self.H, self.K, self.L + ) + + wrong = np.abs(common) ** 2 + for layer, probability in zip( + states.layer_numbers, states.probabilities + ): + amplitude = _poisson_oracle.flat_height_crystal( + int(layer), w_base=w_base + ).F(self.H, self.K, self.L) + wrong = wrong + probability * np.abs(amplitude - common) ** 2 + + correct = model.F2(self.H, self.K, self.L) + self.assertGreater( + np.max(np.abs(correct - wrong)) / np.max(correct), 0.1 + ) + + def test_outer_domain_transforms_reach_each_state(self): + """Every state carries the crystal's area, weight, and domains.""" + profile = self.PROFILES["fractional step"] + transforms = [ + (np.identity(3), 0.65), + (np.diag([1.0, 1.0, 0.5]), 0.35), + ] + reference_uc = _poisson_oracle.layered_cell(2, "reference", a=6.0) + keyargs = {"reference_uc": reference_uc, "domains": transforms} + + model, _, surface, w_base = self.model( + profile, kappa=1.0, weights=[0.4, 0.4], **keyargs + ) + expected = self.oracle_incoherent( + surface, w_base, weights=[0.4], **keyargs + ) + scaled = model.F2(self.H, self.K, self.L) + np.testing.assert_allclose(scaled, expected, rtol=1e-12) + + plain, _, _, _ = self.model(profile, kappa=1.0) + unscaled = plain.F2(self.H, self.K, self.L) + self.assertGreater( + np.max(np.abs(scaled - unscaled)) / np.max(unscaled), 0.1 + ) + + def test_multi_layer_cycles(self): + """The model follows the surface through its layer cycle.""" + profile = self.PROFILES["poisson growth"] + for n_layers in (1, 3): + with self.subTest(n_layers=n_layers): + model, _, surface, w_base = self.model( + profile, kappa=1.0, n_layers=n_layers + ) + np.testing.assert_allclose( + model.F2(self.H, self.K, self.L), + self.oracle_incoherent( + surface, w_base, n_layers=n_layers + ), + rtol=1e-12, + ) + + +class TestPoissonHeightDomainsContract(PoissonModelMixin): + """Construction, validation, and configuration of the Poisson model.""" + + PROFILE = PoissonProfile(1.5, alpha=0.5, tail_probability=1e-14) + + def test_no_F_method_is_defined(self): + """Requesting an amplitude from a mixed state is a category error.""" + model, _, _, _ = self.model(self.PROFILE) + self.assertFalse(hasattr(model, "F")) + + def test_the_fraction_is_a_view_on_the_local_basis(self): + """A fixed and a fitted fraction are the same stored value.""" + model, _, _, _ = self.model(self.PROFILE, kappa=0.35) + self.assertAlmostEqual(model.incoherent_fraction, 0.35) + self.assertAlmostEqual(model.basis[0], 0.35) + + model.incoherent_fraction = 0.6 + self.assertAlmostEqual(model.basis[0], 0.6) + + model.basis[0] = 0.2 + self.assertAlmostEqual(model.incoherent_fraction, 0.2) + + for bad in (-0.1, 1.2, np.nan): + with self.subTest(value=bad): + with self.assertRaisesRegex( + ValueError, "incoherent_fraction" + ): + model.incoherent_fraction = bad + + def test_the_fraction_becomes_a_fit_parameter_on_request(self): + """A setting stays fixed unless it is added as a fit parameter.""" + model, _, _, _ = self.model(self.PROFILE, kappa=0.35) + self.assertEqual(model.n_local_parameters, 0) + + model.addFitParameter( + "incoherent_fraction", + limits=(0.0, 1.0), + name="surface incoherent_fraction", + ) + self.assertEqual(model.fitparnames[0], "surface incoherent_fraction") + np.testing.assert_allclose(model.getInitialParameters()[0], 0.35) + + def test_registered_under_its_stable_key(self): + """The registry resolves a saved type name to this class.""" + models = available_incoherent_models() + self.assertIn("poisson_height_domains", models) + info = models["poisson_height_domains"] + self.assertIs(info.model_class, CTRincoherent.PoissonHeightDomains) + self.assertEqual(info.output_quantity, "F2") + + def test_builtin_registration_cannot_be_replaced(self): + """A built-in key is not available to a third-party model.""" + with self.assertRaisesRegex(ValueError, "built-in"): + + @register_incoherent_model + class _Impostor(CTRincoherent.PoissonHeightDomains): + """Attempts to take over a built-in key.""" + + def test_config_round_trip_keeps_settings_and_fraction(self): + """Settings and the local basis survive a dictionary round-trip.""" + model, crystal, _, _ = self.model(self.PROFILE, kappa=0.35) + model.exact_layer_count = 7 + config = model.to_config() + self.assertEqual(config["type"], "poisson_height_domains") + self.assertEqual( + config["settings"], + {"surface": "surface", "exact_layer_count": 7}, + ) + + restored = create_incoherent_model_from_config(crystal, config) + self.assertEqual(restored.surface, "surface") + self.assertEqual(restored.exact_layer_count, 7) + self.assertAlmostEqual(restored.incoherent_fraction, 0.35) + np.testing.assert_allclose( + restored.F2(self.H, self.K, self.L), + model.F2(self.H, self.K, self.L), + rtol=1e-14, + ) + + def test_rejects_a_missing_or_ambiguous_target(self): + """Selection is by stable name, which must identify one component.""" + w_base = _poisson_oracle.minimum_base_width(self.PROFILE) + crystal, _ = _poisson_oracle.poisson_crystal( + self.PROFILE, w_base=w_base + ) + with self.assertRaisesRegex(ValueError, "No component named"): + CTRincoherent.PoissonHeightDomains(crystal, surface="absent") + + crystal.uc_surface_list[0].name = "surface" + with self.assertRaisesRegex(ValueError, "are named"): + CTRincoherent.PoissonHeightDomains(crystal, surface="surface") + + def test_rejects_a_target_which_is_not_a_poisson_surface(self): + """The model only knows how to enumerate Poisson heights.""" + w_base = _poisson_oracle.minimum_base_width(self.PROFILE) + crystal, _ = _poisson_oracle.poisson_crystal( + self.PROFILE, w_base=w_base + ) + with self.assertRaisesRegex(ValueError, "not a PoissonSurface"): + CTRincoherent.PoissonHeightDomains(crystal, surface="film") + + def test_rejects_a_component_stacked_above_the_target(self): + """A component above the surface may move with the chosen height.""" + w_base = _poisson_oracle.minimum_base_width(self.PROFILE) + crystal, surface = _poisson_oracle.poisson_crystal( + self.PROFILE, w_base=w_base + ) + cap = CTRcalc.UnitCell([3.0, 3.0, 6.0], [90.0, 90.0, 90.0], name="cap") + cap.addAtom("O", [0.0, 0.0, 0.0], 0.1, 0.1, 1.0, layer=0) + stacked = CTRcalc.SXRDCrystal( + _poisson_oracle.layered_cell(2, "bulk"), + crystal.uc_surface_list[0], + surface, + cap, + stacking=np.array([1, 2, 3]), + ) + with self.assertRaisesRegex(ValueError, "topmost component"): + CTRincoherent.PoissonHeightDomains(stacked, surface="surface") + + def test_rejects_an_invalid_state_count(self): + """``exact_layer_count`` is a positive integer policy setting.""" + w_base = _poisson_oracle.minimum_base_width(self.PROFILE) + crystal, _ = _poisson_oracle.poisson_crystal( + self.PROFILE, w_base=w_base + ) + for bad in (0, -3, 1.5): + with self.subTest(exact_layer_count=bad): + with self.assertRaisesRegex(ValueError, "positive integer"): + CTRincoherent.PoissonHeightDomains( + crystal, surface="surface", exact_layer_count=bad + ) + + def test_dwba_is_refused_before_any_prediction(self): + """An F2 wrapper must never be reinterpreted as a reflectivity.""" + model, _, _, _ = self.model(self.PROFILE) + model.validate(forward_model="kinematical") + with self.assertRaisesRegex(ValueError, "does not support"): + model.validate(forward_model="dwba") + + def test_survives_the_deep_copy_the_optimizer_performs(self): + """Component selection by name must outlive a copy.""" + model, _, _, _ = self.model(self.PROFILE, kappa=0.4) + copied = copy.deepcopy(model) + self.assertIsNot(copied.coherent_model, model.coherent_model) + np.testing.assert_allclose( + copied.F2(self.H, self.K, self.L), + model.F2(self.H, self.K, self.L), + rtol=1e-14, + ) + + if __name__ == "__main__": unittest.main() From ccd56df93fb86747296a2bd082ef509eacd793dc Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Fri, 11 Sep 2026 20:59:15 -0400 Subject: [PATCH 11/13] feat(phys): fit incoherent models through CTROptimizer Resolution now acts on F2 before the square root, and an incoherent wrapper is accepted through the existing model argument. optimizer.xtal stays the coherent crystal, so callbacks and constraints are unchanged. --- orgui/datautils/xrayutils/CTRopt.py | 136 +++++- orgui/datautils/xrayutils/CTRresolution.py | 40 +- orgui/datautils/xrayutils/test/test_CTRopt.py | 401 +++++++++++++++++- .../xrayutils/test/test_CTRresolution.py | 67 +++ 4 files changed, 615 insertions(+), 29 deletions(-) diff --git a/orgui/datautils/xrayutils/CTRopt.py b/orgui/datautils/xrayutils/CTRopt.py index 4e7a2d7..bbc0c23 100644 --- a/orgui/datautils/xrayutils/CTRopt.py +++ b/orgui/datautils/xrayutils/CTRopt.py @@ -40,6 +40,7 @@ from .. import util from .CTRcalc import SXRDCrystal +from .CTRincoherent import IncoherentModel from . import CTRplotutil, CTRresolution @dataclass(frozen=True) @@ -86,9 +87,30 @@ class _ScaleEstimationError(ValueError): class CTROptimizer: def __init__(self, xtal, CTRs, *, scale_policy=None): + """Fit one forward model against a CTR collection. + + :param xtal: + Forward model. Either a coherent ``SXRDCrystal`` or an + ``IncoherentF2Model`` wrapping one. The model is deep-copied once, + so later changes to the caller's object do not reach the fit. + :param CTRplotutil.CTRCollection CTRs: + Measured CTRs. + :param scale_policy: + Optional initial scale policy mapping. + """ self.CTRs = copy.deepcopy(CTRs) self.CTRs.sort(key=lambda x: abs(x.hk[0]) + abs(x.hk[1])) - self.xtal = copy.deepcopy(xtal) + # `self.model` is the fitted forward model and `self.xtal` stays the + # primary coherent crystal, so registered callbacks, displacement + # constraints, and public `optimizer.xtal` access keep receiving an + # SXRDCrystal whether or not the fit is wrapped. + self.model = copy.deepcopy(xtal) + if isinstance(self.model, IncoherentModel): + self.xtal = self.model.coherent_model + else: + self.xtal = self.model + self.n_parameters = None + self._prepared_signature = None self._scale_policy_defaults = { "structure_factor": "scaled", "reflectivity": "fixed", @@ -357,6 +379,28 @@ def _validate_kinematical_input(self): f"{ctr!r}: kinematical CTR fitting supports " "structure-factor data only." ) + if isinstance(self.model, IncoherentModel): + if self.model.output_quantity != "F2": + raise ValueError( + f"{self.model.model_type} returns " + f"{self.model.output_quantity!r}, which the kinematical " + "path cannot consume as a squared structure factor" + ) + self.model.validate(forward_model="kinematical") + + def _kinematic_F2(self, h, k, l): # noqa: N802,E741 + """Return the forward model's squared structure factor. + + Prefers ``F2``. Crystal-like objects and test doubles which implement + only ``F`` keep working through the squared-modulus fallback. + + :returns: + Real squared structure factor in the squared units of ``F``. + """ + squared = getattr(self.model, "F2", None) + if callable(squared): + return np.asarray(squared(h, k, l), dtype=np.float64) + return np.abs(self.model.F(h, k, l)) ** 2 def _validate_measurement_input(self): """Require aligned finite real data and positive uncertainties.""" @@ -392,6 +436,16 @@ def _validate_measurement_input(self): def _validate_dwba_input(self): """Validate measurement metadata required by live DWBA prediction.""" + if isinstance(self.model, IncoherentModel): + # Rejected before any DWBA code runs. Changing the flat surface + # height changes the optical reference profile and its internal + # fields, so a kinematical F2 mixture is not a reflectivity and + # must never be reinterpreted as one. + raise ValueError( + f"{self.model.model_type} does not support the DWBA forward " + f"model: it returns {self.model.output_quantity}, not a " + "reflectivity" + ) if not isinstance(self.xtal, SXRDCrystal): raise TypeError("DWBA fitting requires an SXRDCrystal model") if self._fit_resolution and self.resolution_calculation == "sample": @@ -578,13 +632,18 @@ def _update_resolution_cache(self): return if self.resolution_calculation == "sample": self._resolution_calculated_ctrs = CTRresolution.sample_structure_factor( - self.CTRs, self.xtal, self.resolution + self.CTRs, self.model, self.resolution ) return input_ctrs = self._resolution_input_collection() + # `fast_convolve` squares its input, convolves, and takes one square + # root, so feeding it sqrt(F2) applies the resolution to F2 before the + # square root without a second collection-building path. for source, calculated in zip(self.CTRs, input_ctrs): - calculated.sfI = np.abs(self.xtal.F(source.harr, source.karr, source.l)) + calculated.sfI = np.sqrt( + self._kinematic_F2(source.harr, source.karr, source.l) + ) self._resolution_calculated_ctrs = CTRresolution.fast_convolve( input_ctrs, self.resolution ) @@ -810,7 +869,7 @@ def _calculated_value(self, ctr, index): if self._dwba_enabled: return self._dwba_prediction(ctr) if self.resolution is None: - return np.abs(self.xtal.F(ctr.harr, ctr.karr, ctr.l)) + return np.sqrt(self._kinematic_F2(ctr.harr, ctr.karr, ctr.l)) if self._resolution_calculated_ctrs is None: self._update_resolution_cache() return self._resolution_calculated_ctrs[index].sfI @@ -973,24 +1032,28 @@ def _evaluate(self, x=None): return calculations def _model_parameters(self): - """Return the crystal-owned tail of the fit parameter vector.""" - return self.xtal.getInitialParameters() + """Return the model-owned tail of the fit parameter vector. + + For a wrapper this is its local block followed by the wrapped + crystal's, which the wrapper splits internally. + """ + return self.model.getInitialParameters() def _set_model_parameters(self, parameters): - """Set parameters owned by the fitted crystal model.""" - self.xtal.setParameters(parameters) + """Set parameters owned by the fitted forward model.""" + self.model.setParameters(parameters) def _prepend_model_bounds(self, bounds): - """Add subclass-owned parameters ahead of the crystal bounds.""" + """Add subclass-owned parameters ahead of the model bounds.""" return bounds def _set_model_errors(self, errors): - """Forward the model-owned error tail to the crystal.""" - self.xtal.setFitErrors(errors) + """Forward the model-owned error tail to the forward model.""" + self.model.setFitErrors(errors) def _model_parameter_names(self): - """Return names for subclass and crystal model parameters.""" - return list(self.xtal.fitparnames) + """Return names for subclass and forward-model parameters.""" + return list(self.model.fitparnames) def _fit_parameter_names(self): """Build names in the same order as the fit parameter vector.""" @@ -1037,8 +1100,11 @@ def prepareFit(self): self._validate_dwba_input() else: self._validate_kinematical_input() + # `startp`, `lower_bounds` and `higher_bounds` keep their existing + # model-block scope. `self.bounds`, `fitparnames`, + # `get_parameters()` and `n_parameters` describe the full vector. self.startp, self.lower_bounds, self.higher_bounds = ( - self.xtal.getStartParamAndLimits() + self.model.getStartParamAndLimits() ) self.bounds = self._prepend_model_bounds( (self.lower_bounds, self.higher_bounds) @@ -1056,14 +1122,53 @@ def prepareFit(self): else: self.nic = 0 self.fitparnames = self._fit_parameter_names() - self.priors = self.xtal.priors + # Model-scoped, matching current behavior: resolution, callback + # and angle-correction prefixes expose no prior API. + self.priors = self.model.priors self._prepared = True + self.n_parameters = len(self.get_parameters()) + if len(self.fitparnames) != self.n_parameters or len( + self.bounds[0] + ) != self.n_parameters: + raise ValueError( + "Prepared fit vector is inconsistent: " + f"{self.n_parameters} values, {len(self.fitparnames)} " + f"names, {len(self.bounds[0])} bounds" + ) + self._prepared_signature = self._model_signature() self._evaluate() except Exception: self._prepared = False self._invalidate_calculated_results() raise + def _model_signature(self): + """Return the structural layout of the model parameter block. + + Names, count, and bounds. A change to any of them means the prepared + vector no longer describes the model, so slicing it would silently + address the wrong parameters. + """ + start, lower, upper = self.model.getStartParamAndLimits() + return ( + tuple(self.model.fitparnames), + len(start), + tuple(np.asarray(lower, dtype=np.float64).tolist()), + tuple(np.asarray(upper, dtype=np.float64).tolist()), + ) + + def _require_current_signature(self): + """Reject a stale layout after structural parameter changes.""" + if self._prepared_signature is None: + return + if self._model_signature() != self._prepared_signature: + self._prepared = False + raise ValueError( + "Fit parameters were added, removed, renamed, or re-bounded " + "on optimizer.model after prepareFit(). Call prepareFit() " + "again before evaluating or setting parameters." + ) + def get_bounds(self): return self.bounds @@ -1079,6 +1184,7 @@ def get_parameters(self): def set_parameters(self, x): """Set resolution, callback, and model parameters in layout order.""" + self._require_current_signature() self._invalidate_calculated_results() x = self._split_resolution_parameters(x) counter = 0 diff --git a/orgui/datautils/xrayutils/CTRresolution.py b/orgui/datautils/xrayutils/CTRresolution.py index 1b8029f..4616642 100644 --- a/orgui/datautils/xrayutils/CTRresolution.py +++ b/orgui/datautils/xrayutils/CTRresolution.py @@ -497,7 +497,9 @@ def sample_structure_factor(ctrs, crystal, resolution, quadrature_order=25): :param CTRCollection ctrs: Collection supplying the requested HKL points and optional angles. :param crystal: - Crystal-like object providing ``F(h, k, l)``. + Forward model providing ``F2(h, k, l)``, or a crystal-like object + providing only ``F(h, k, l)``. ``F2`` is preferred, because a mixed + state has a squared structure factor but no unique complex amplitude. :param ResolutionFunction resolution: Box or Gaussian L-resolution function. :param int quadrature_order: @@ -505,15 +507,20 @@ def sample_structure_factor(ctrs, crystal, resolution, quadrature_order=25): 25. :returns: A new collection containing effective amplitudes - ``sqrt(integrated(abs(F)**2))``. + ``sqrt(integrated(F2))``. :rtype: CTRCollection """ if not isinstance(ctrs, CTRCollection): raise TypeError("ctrs must be a CTRCollection") if not isinstance(resolution, ResolutionFunction): raise TypeError("resolution must be a ResolutionFunction") - if not hasattr(crystal, "F") or not callable(crystal.F): - raise TypeError("crystal must provide a callable F(h, k, l) method") + squared = getattr(crystal, "F2", None) + if not callable(squared) and not ( + hasattr(crystal, "F") and callable(crystal.F) + ): + raise TypeError( + "crystal must provide a callable F2(h, k, l) or F(h, k, l) method" + ) quadrature_order = _validate_quadrature_order(quadrature_order) sampled = [] @@ -524,21 +531,28 @@ def sample_structure_factor(ctrs, crystal, resolution, quadrature_order=25): angles = getattr(ctr, "angles", None) def structure_factor_intensity(h_samples, k_samples, l_samples): - structure_factor = np.asarray( - crystal.F(h_samples, k_samples, l_samples) - ) - try: - structure_factor = np.broadcast_to( - structure_factor, h_samples.shape + if callable(squared): + values = np.asarray(squared(h_samples, k_samples, l_samples)) + source = "crystal.F2" + else: + values = ( + np.abs( + np.asarray( + crystal.F(h_samples, k_samples, l_samples) + ) + ) + ** 2 ) + source = "crystal.F" + try: + intensity = np.broadcast_to(values, h_samples.shape) except ValueError as exc: raise ValueError( - "crystal.F returned an incompatible array shape" + f"{source} returned an incompatible array shape" ) from exc - intensity = np.abs(structure_factor) ** 2 if not np.all(np.isfinite(intensity)): raise ValueError( - "crystal.F returned non-finite structure factors" + f"{source} returned non-finite structure factors" ) return intensity diff --git a/orgui/datautils/xrayutils/test/test_CTRopt.py b/orgui/datautils/xrayutils/test/test_CTRopt.py index cfc6fd8..0aa6d5f 100644 --- a/orgui/datautils/xrayutils/test/test_CTRopt.py +++ b/orgui/datautils/xrayutils/test/test_CTRopt.py @@ -14,7 +14,16 @@ import numpy as np import pytest -from .. import CTRcalc, CTRopt, CTRplotutil, CTRresolution, CTRsymmetry +from .. import ( + CTRcalc, + CTRdistributions, + CTRincoherent, + CTRopt, + CTRplotutil, + CTRresolution, + CTRsymmetry, +) +from . import _poisson_oracle def _angles(gamma, delta=None): @@ -1708,6 +1717,396 @@ def test_unprefixed_layout_is_unchanged(self): self.assertEqual(len(optimizer.xtal.error_calls), 1) +class _F2OnlyModel: + """Forward model exposing only ``F2``, like an incoherent wrapper. + + Kept distinct from ``FitCrystal`` so the optimizer's preference for + ``F2`` is visible: this deliberately has no ``F`` at all. + """ + + def __init__(self, parameters=(1.0,)): + self.parameters = np.asarray(parameters, dtype=np.float64) + self.fitparnames = [f"xtal_{i}" for i in range(self.parameters.size)] + self.priors = [] + self.errors = None + self.error_calls = [] + self.offset = 1.0 + self.slope = 0.5 + + def F2(self, h, k, l): # noqa: N802,E741 + """Return the squared structure factor directly.""" + lvalues = np.asarray(l, dtype=np.float64) + value = self.offset + self.slope * lvalues + for i, parameter in enumerate(self.parameters): + value = value + parameter * lvalues ** (i + 2) + return value**2 + + def getStartParamAndLimits(self): # noqa: N802 + """Return the start parameters and their lower and upper limits.""" + return ( + self.parameters.copy(), + np.full(self.parameters.size, 0.1), + np.full(self.parameters.size, 10.0), + ) + + def getInitialParameters(self): # noqa: N802 + """Return the current fit parameters.""" + return self.parameters.copy() + + def setParameters(self, parameters): # noqa: N802 + """Set the fit parameters.""" + self.parameters = np.asarray(parameters, dtype=np.float64) + + def setFitErrors(self, errors): # noqa: N802 + """Store the parameter errors and record the call.""" + self.errors = None if errors is None else np.asarray(errors) + self.error_calls.append(self.errors) + + +class TestForwardModelQuantity(unittest.TestCase): + """The optimizer consumes ``F2`` and falls back to ``abs(F) ** 2``.""" + + def test_an_F_only_model_still_works(self): + """Existing crystal-like objects implement only ``F``.""" + optimizer = CTRopt.CTROptimizer(FitCrystal(), _fixture_ctrs()) + optimizer.prepareFit() + crystal = optimizer.model + expected = np.abs( + crystal.F( + optimizer.CTRs[0].harr, + optimizer.CTRs[0].karr, + optimizer.CTRs[0].l, + ) + ) + np.testing.assert_allclose( + optimizer._calculated_value(optimizer.CTRs[0], 0), expected + ) + + def test_an_F2_model_is_preferred(self): + """A model exposing ``F2`` is used without ever calling ``F``.""" + optimizer = CTRopt.CTROptimizer(_F2OnlyModel(), _fixture_ctrs()) + optimizer.prepareFit() + ctr = optimizer.CTRs[0] + np.testing.assert_allclose( + optimizer._calculated_value(ctr, 0), + np.sqrt(optimizer.model.F2(ctr.harr, ctr.karr, ctr.l)), + ) + + def test_the_model_and_the_crystal_coincide_without_a_wrapper(self): + """``optimizer.xtal`` keeps its meaning for a coherent fit.""" + optimizer = CTRopt.CTROptimizer(FitCrystal(), _fixture_ctrs()) + self.assertIs(optimizer.model, optimizer.xtal) + + +class TestPreparedVectorLayout(unittest.TestCase): + """``n_parameters`` describes the full vector; ``startp`` does not.""" + + @staticmethod + def _optimizer(callbacks, fit_resolution, angle_correction): + """Build one optimizer with the requested optional prefixes.""" + factory = ( + CTRopt.CTROptAngleCorrection + if angle_correction + else CTRopt.CTROptimizer + ) + optimizer = factory(FitCrystal((1.0, 2.0)), _fixture_ctrs()) + if angle_correction: + optimizer.useAnglecorr = True + if fit_resolution: + optimizer.fit_resolution( + CTRresolution.BoxResolution(0.1, 0.1, 0.1), + lower_bounds=[0.0, 0.0, 0.0], + higher_bounds=[1.0, 1.0, 1.0], + ) + for index in range(callbacks): + _register_callback(optimizer, f"cb_{index}", 1) + return optimizer + + def test_full_vector_agrees_with_n_parameters(self): + """Values, names, and bounds all describe the same prepared vector.""" + for callbacks in (0, 2): + for fit_resolution in (False, True): + for angle_correction in (False, True): + with self.subTest( + callbacks=callbacks, + fit_resolution=fit_resolution, + angle_correction=angle_correction, + ): + optimizer = self._optimizer( + callbacks, fit_resolution, angle_correction + ) + optimizer.prepareFit() + total = optimizer.n_parameters + self.assertEqual(len(optimizer.get_parameters()), total) + self.assertEqual(len(optimizer.fitparnames), total) + self.assertEqual(len(optimizer.get_bounds()[0]), total) + self.assertEqual(len(optimizer.get_bounds()[1]), total) + + expected = ( + 2 + + callbacks + + (3 if fit_resolution else 0) + + (2 if angle_correction else 0) + ) + self.assertEqual(total, expected) + + def test_startp_keeps_its_model_block_scope(self): + """``startp`` is the preparation snapshot of the model block only. + + Redefining it as the full vector would silently change its length for + every existing fit using callbacks or fitted resolution. + """ + optimizer = self._optimizer(2, True, True) + optimizer.prepareFit() + self.assertEqual( + len(optimizer.startp), + len(optimizer.model.getInitialParameters()), + ) + self.assertLess(len(optimizer.startp), optimizer.n_parameters) + self.assertEqual(len(optimizer.lower_bounds), len(optimizer.startp)) + self.assertEqual(len(optimizer.higher_bounds), len(optimizer.startp)) + + def test_get_parameters_is_the_live_vector(self): + """Later value changes appear without re-preparing.""" + optimizer = self._optimizer(0, False, False) + optimizer.prepareFit() + snapshot = optimizer.startp.copy() + optimizer.set_parameters([3.0, 4.0]) + np.testing.assert_allclose(optimizer.get_parameters(), [3.0, 4.0]) + np.testing.assert_allclose(optimizer.startp, snapshot) + + def test_structural_changes_after_preparation_are_rejected(self): + """A stale layout must not be sliced silently.""" + optimizer = CTRopt.CTROptimizer(FitCrystal((1.0, 2.0)), _fixture_ctrs()) + optimizer.prepareFit() + optimizer.model.fitparnames = ["xtal_0", "renamed"] + with self.assertRaisesRegex(ValueError, "prepareFit"): + optimizer.set_parameters([1.0, 2.0]) + + +class TestIncoherentWrapperIntegration(unittest.TestCase): + """An incoherent wrapper fits through the existing model argument.""" + + def build(self, kappa=0.35, fit_fraction=True): + """Return a wrapped model and its coherent crystal.""" + profile = CTRdistributions.PoissonProfile( + 1.5, alpha=0.5, tail_probability=1e-12 + ) + w_base = _poisson_oracle.minimum_base_width(profile) + crystal, _ = _poisson_oracle.poisson_crystal(profile, w_base=w_base) + crystal.apply_stacking() + crystal["film"].addFitParameter("W", limits=(1.0, 10.0)) + model = CTRincoherent.PoissonHeightDomains( + crystal, surface="surface", incoherent_fraction=kappa + ) + if fit_fraction: + model.addFitParameter( + "incoherent_fraction", + limits=(0.0, 1.0), + name="surface incoherent_fraction", + ) + return model + + def optimizer(self, model, **keyargs): + """Return a prepared optimizer over the shared CTR fixture.""" + fit = CTRopt.CTROptimizer( + model, _fixture_ctrs(), scale_policy={"F": "fixed"}, **keyargs + ) + return fit + + def test_xtal_stays_the_coherent_crystal(self): + """Callbacks and constraints must keep receiving an SXRDCrystal.""" + fit = self.optimizer(self.build()) + self.assertIsInstance(fit.model, CTRincoherent.PoissonHeightDomains) + self.assertIs(fit.xtal, fit.model.coherent_model) + self.assertIsInstance(fit.xtal, CTRcalc.SXRDCrystal) + + def test_the_model_is_deep_copied(self): + """Mutating the caller's model must not reach the fit.""" + model = self.build() + fit = self.optimizer(model) + self.assertIsNot(fit.model, model) + model.incoherent_fraction = 0.9 + self.assertAlmostEqual(fit.model.incoherent_fraction, 0.35) + + def test_parameter_layout_is_local_then_crystal(self): + """The wrapper's block sits ahead of the crystal's, after prefixes.""" + fit = self.optimizer(self.build()) + _register_callback(fit, "cb", 1) + fit.fit_resolution( + CTRresolution.BoxResolution(0.1, 0.1, 0.1), + lower_bounds=[0.0, 0.0, 0.0], + higher_bounds=[1.0, 1.0, 1.0], + ) + fit.prepareFit() + + names = list(fit.fitparnames) + self.assertEqual(names[:3], [ + "resolution_delta_l_0", + "resolution_delta_l_1", + "resolution_delta_l_2", + ]) + self.assertEqual(names[4], "surface incoherent_fraction") + self.assertEqual(names[5:], list(fit.xtal.fitparnames)) + self.assertEqual(fit.n_parameters, len(names)) + + def test_priors_stay_model_scoped(self): + """Priors cover the wrapper's local and coherent block only.""" + model = self.build() + fit = self.optimizer(model) + fit.fit_resolution( + CTRresolution.BoxResolution(0.1, 0.1, 0.1), + lower_bounds=[0.0, 0.0, 0.0], + higher_bounds=[1.0, 1.0, 1.0], + ) + fit.prepareFit() + self.assertEqual(len(fit.priors), len(fit.model.priors)) + self.assertLess(len(fit.priors), fit.n_parameters) + + def test_direct_prediction_is_the_square_root_of_F2(self): + """Without resolution the stored prediction is ``sqrt(F2)``.""" + fit = self.optimizer(self.build()) + fit.prepareFit() + for index, ctr in enumerate(fit.CTRs): + np.testing.assert_allclose( + fit._calculated_value(ctr, index), + np.sqrt(fit.model.F2(ctr.harr, ctr.karr, ctr.l)), + rtol=1e-12, + ) + + def test_sampled_resolution_integrates_F2(self): + """Quadrature evaluates ``F2`` at every sampled coordinate.""" + fit = self.optimizer(self.build()) + resolution = CTRresolution.BoxResolution(0.05, 0.0, 0.0) + fit.set_resolution(resolution, calculation="sample") + fit.prepareFit() + + expected = CTRresolution.sample_structure_factor( + fit.CTRs, fit.model, resolution + ) + for index, ctr in enumerate(fit.CTRs): + np.testing.assert_allclose( + fit._calculated_value(ctr, index), expected[index].sfI, + rtol=1e-12, + ) + + def test_fast_convolution_broadens_F2_before_the_square_root(self): + """The broadened prediction is ``sqrt(convolved(F2))``. + + The dense rod matters: with only a few L points the kernel has almost + nothing to act on, and broadening before or after the square root + would be indistinguishable. + """ + lvalues = np.linspace(0.3, 2.7, 25) + dense = CTRplotutil.CTRCollection( + [ + CTRplotutil.CTR( + (0.0, 0.0), + lvalues.copy(), + np.ones_like(lvalues), + np.full_like(lvalues, 0.1), + ) + ] + ) + fit = CTRopt.CTROptimizer( + self.build(), dense, scale_policy={"F": "fixed"} + ) + resolution = CTRresolution.BoxResolution(0.3, 0.0, 0.0) + fit.set_resolution(resolution, calculation="convolve") + fit.prepareFit() + + ctr = fit.CTRs[0] + squared = fit.model.F2(ctr.harr, ctr.karr, ctr.l) + correct = np.sqrt( + CTRresolution.fast_convolve_intensity( + ctr.harr, ctr.karr, ctr.l, squared, resolution + ) + ) + # The wrong order: broaden the amplitude and square afterwards. + wrong = CTRresolution.fast_convolve_intensity( + ctr.harr, ctr.karr, ctr.l, np.sqrt(squared), resolution + ) + + np.testing.assert_allclose( + fit._calculated_value(ctr, 0), correct, rtol=1e-12 + ) + self.assertGreater( + np.max(np.abs(correct - wrong)) / np.max(correct), 1e-4 + ) + self.assertGreater( + np.max(np.abs(correct - np.sqrt(squared))) / np.max(correct), 1e-4 + ) + + def test_predictions_track_the_mixing_fraction(self): + """A coherent wrapper predicts what the bare crystal predicts.""" + model = self.build(kappa=0.0) + fit = self.optimizer(model) + fit.prepareFit() + coherent = np.copy(fit.flat_prediction()) + + # Same crystal, no wrapper, same scale policy: identical predictions. + bare = CTRopt.CTROptimizer( + model.coherent_model, + _fixture_ctrs(), + scale_policy={"F": "fixed"}, + ) + bare.prepareFit() + np.testing.assert_allclose( + coherent, bare.flat_prediction(), rtol=1e-10 + ) + + values = fit.get_parameters() + values[0] = 1.0 + fit.set_parameters(values) + incoherent = np.copy(fit.flat_prediction()) + self.assertGreater( + np.max(np.abs(incoherent - coherent)) / np.max(coherent), 1e-3 + ) + + def test_error_slices_reach_their_owners(self): + """Resolution, callback, wrapper-local, and crystal errors split.""" + fit = self.optimizer(self.build()) + callback, _ = _register_callback(fit, "cb", 2) + fit.fit_resolution( + CTRresolution.BoxResolution(0.1, 0.1, 0.1), + lower_bounds=[0.0, 0.0, 0.0], + higher_bounds=[1.0, 1.0, 1.0], + ) + fit.prepareFit() + + errors = 0.01 * (1.0 + np.arange(fit.n_parameters)) + fit.set_errors(errors) + + np.testing.assert_allclose(fit.resolution_errors, errors[:3]) + np.testing.assert_allclose(callback.errors, errors[3:5]) + np.testing.assert_allclose(fit.model.getFitErrors(), errors[5:]) + + def test_callbacks_receive_the_crystal(self): + """A callback signature takes an SXRDCrystal, not the wrapper.""" + fit = self.optimizer(self.build()) + seen = [] + + def apply(xtal, x): + seen.append(xtal) + + fit.register_fit_callback(apply, [0.1], [5.0], [1.0], name="cb") + fit.prepareFit() + # Callbacks run when parameters are applied, not during preparation. + fit.set_parameters(fit.get_parameters()) + self.assertTrue(seen) + for received in seen: + self.assertIs(received, fit.xtal) + + def test_dwba_is_rejected_during_preparation(self): + """An F2 wrapper must fail before any DWBA code runs.""" + fit = self.optimizer(self.build()) + fit.set_dwba(True) + with self.assertRaisesRegex( + ValueError, "does not support the DWBA forward model" + ): + fit.prepareFit() + + class TestCallbackBounds(unittest.TestCase): """Increment 1: the ``FitCallback`` bounds ordering repair.""" diff --git a/orgui/datautils/xrayutils/test/test_CTRresolution.py b/orgui/datautils/xrayutils/test/test_CTRresolution.py index 36502cb..972828c 100644 --- a/orgui/datautils/xrayutils/test/test_CTRresolution.py +++ b/orgui/datautils/xrayutils/test/test_CTRresolution.py @@ -573,5 +573,72 @@ def test_collection_forwards_all_z_mode_constraints(self): np.testing.assert_allclose(ctrs[1].angles["gamma"], [0.25]) +class _AmplitudeOnly: + """Crystal-like object exposing only ``F``.""" + + def F(self, h, k, l): # noqa: N802,E741 + """Return a simple complex amplitude.""" + return np.asarray(l, dtype=np.float64) + 2.0j + + +class _SquaredAndAmplitude(_AmplitudeOnly): + """Model whose ``F2`` deliberately disagrees with ``abs(F) ** 2``. + + Only a model which prefers ``F2`` can tell the two apart, which is what + makes the preference observable rather than assumed. + """ + + def F2(self, h, k, l): # noqa: N802,E741 + """Return a squared structure factor unrelated to ``abs(F) ** 2``.""" + return 7.0 * np.ones_like(np.asarray(l, dtype=np.float64)) + + +class TestSampleStructureFactorQuantity(unittest.TestCase): + """``sample_structure_factor`` prefers ``F2`` over ``abs(F) ** 2``.""" + + def collection(self): + """Return a single rod with room for the kernel to act.""" + lvalues = np.linspace(0.5, 2.5, 9) + ctr = CTRplotutil.CTR( + (0.0, 0.0), + lvalues, + np.ones_like(lvalues), + np.full_like(lvalues, 0.1), + ) + return CTRplotutil.CTRCollection([ctr]) + + def test_F2_is_used_when_available(self): + """A constant ``F2`` broadens to its own square root.""" + ctrs = self.collection() + resolution = CTRresolution.BoxResolution(0.2, 0.0, 0.0) + result = CTRresolution.sample_structure_factor( + ctrs, _SquaredAndAmplitude(), resolution + ) + np.testing.assert_allclose( + result[0].sfI, np.sqrt(7.0), rtol=1e-9 + ) + + def test_amplitude_only_models_still_work(self): + """Objects implementing only ``F`` keep the squared-modulus path.""" + ctrs = self.collection() + resolution = CTRresolution.BoxResolution(0.2, 0.0, 0.0) + amplitude = CTRresolution.sample_structure_factor( + ctrs, _AmplitudeOnly(), resolution + ) + self.assertTrue(np.all(np.isfinite(amplitude[0].sfI))) + self.assertFalse( + np.allclose(amplitude[0].sfI, np.sqrt(7.0), rtol=1e-9) + ) + + def test_a_model_with_neither_method_is_rejected(self): + """The quantity boundary must be explicit.""" + with self.assertRaisesRegex(TypeError, "F2.*or F"): + CTRresolution.sample_structure_factor( + self.collection(), + object(), + CTRresolution.BoxResolution(0.2, 0.0, 0.0), + ) + + if __name__ == "__main__": unittest.main() From d0a19c68e1df39299b4b85db9dad17d8e84556d8 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Sat, 12 Sep 2026 00:59:04 -0400 Subject: [PATCH 12/13] docs(phys): document the incoherent CTR quantities Separates F, F2, r, R and detector counts along the CTR path. Config restore now names the settings which failed to construct the model. --- CHANGELOG.md | 21 +++ doc/source/ctr_structure_factors.rst | 156 +++++++++++++++++- orgui/datautils/xrayutils/CTRincoherent.py | 49 +++++- .../xrayutils/test/test_CTRincoherent.py | 111 +++++++++++++ 4 files changed, 332 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d239056..506db3a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,27 @@ This is the changelog for the software orGUI, written by Timo Fuchs Scientific and analysis additions: +- **Added opt-in incoherent CTR models for large surface height domains.** + A `PoissonSurface` distributes surface heights, and by default those heights + add as amplitudes: the coherent limit, in which every height lies inside one + coherence patch. When the lateral height domains are large compared with the + projected coherence area, each patch instead sees one flat height and the + patches add in intensity. The new `CTRincoherent.PoissonHeightDomains` + interpolates between the two through a dimensionless `incoherent_fraction`, + mixing the *complete* crystal amplitude of each height so that bulk-surface + and Film-surface interference stays inside every domain. Constructing a + `PoissonSurface` does not opt into the averaging: wrapping the crystal is + the explicit opt-in, and existing scripts, saved crystals, and optimizer + setups stay coherent and numerically unchanged. The fraction is a fixed + setting until it is added as a fit parameter, after which the wrapper is + fitted through `CTROptimizer`'s existing model argument, with + `optimizer.xtal` still the coherent crystal for callbacks and constraints. + `SXRDCrystal` gains `F2` and `evaluate_kinematic`, resolution now acts on + `F2` before the conversion back to a stored amplitude, and `CTROptimizer` + gains `n_parameters` for the full prepared vector. The first implementation + is kinematical: combining it with DWBA, or targeting a surface which is not + the topmost component, fails explicitly during `prepareFit`. + - **Added live DWBA predictions to CTR fitting.** ``CTROptimizer.set_dwba`` now evaluates the optimizer-owned crystal through the semi-infinite DWBA model, forms predictions independently in each dataset's stored F or diff --git a/doc/source/ctr_structure_factors.rst b/doc/source/ctr_structure_factors.rst index dbb9c47..2c38a84 100644 --- a/doc/source/ctr_structure_factors.rst +++ b/doc/source/ctr_structure_factors.rst @@ -684,6 +684,132 @@ No illuminated footprint, detector response, or experimental scale factor is included. Calculated intensity is proportional to :math:`|F_{\mathrm{crystal}}|^2`. +Amplitudes, squared structure factors, and reflectivity +-------------------------------------------------------- + +Four quantities appear along the CTR path and are deliberately not +interchangeable. + +``F`` + The coherent complex structure factor returned by + :meth:`~orgui.datautils.xrayutils.CTRcalc.SXRDCrystal.F`, in electrons per + lateral cell of the reference unit cell. It carries a phase. + +``F2`` + The squared structure factor + :meth:`~orgui.datautils.xrayutils.CTRcalc.SXRDCrystal.F2`, equal to + ``abs(F) ** 2`` for a coherent crystal, in the squared units of ``F``. It + is real and nonnegative and carries no phase. + +``r`` + The optical reflection amplitude of the DWBA path, a dimensionless complex + ratio formed from the electron-density profile. + +``R`` + Reflectivity, the dimensionless intensity ratio ``abs(r) ** 2``. + +None of these is detector counts. Incident flux, illuminated footprint, +polarization and Lorentz factors, detector response, acquisition time, +background, and the fitted experimental scale all sit outside them. + +An incoherent model has a ``F2`` but no ``F``: a mixed state has no unique +complex amplitude, so requesting one is a category error rather than a +missing feature. ``F2`` never changes meaning with the optimizer mode; a +reflectivity is returned from a separately named boundary, never from a +method called ``F2``. + +Incoherent height domains +--------------------------------------------- + +A :class:`~orgui.datautils.xrayutils.CTRfilm.PoissonSurface` describes a +distribution of surface heights. By default those heights are added as +*amplitudes*, which is the coherent limit: every height lies inside one +coherence patch, and the height distribution enters through the squared +modulus of its characteristic function. + +When the lateral height domains are large compared with the projected +coherence area, each patch instead sees a single flat height and the patches +add in intensity. The two limits are + +.. math:: + + F^2_{\mathrm{coherent}} = \left| \sum_n p_n A_n \right|^2, + \qquad + F^2_{\mathrm{incoherent}} = \sum_n p_n \left| A_n \right|^2, + +where :math:`p_n` is the probability of height state :math:`n` and +:math:`A_n` is the **complete** crystal amplitude when that height covers the +patch: bulk, Film, the flat-height Film correction, and the exposed +termination together. Averaging only the surface correction would drop the +bulk-surface and Film-surface interference inside each domain and is a +different quantity. + +:class:`~orgui.datautils.xrayutils.CTRincoherent.PoissonHeightDomains` +interpolates between them with a dimensionless ``incoherent_fraction`` +:math:`\kappa` in ``[0, 1]``: + +.. math:: + + F^2_\kappa = (1 - \kappa) F^2_{\mathrm{coherent}} + + \kappa F^2_{\mathrm{incoherent}}. + +For a patch containing :math:`N_{\mathrm{eff}}` independent equal-area +domains, :math:`\kappa \approx 1 / N_{\mathrm{eff}}`. A CTR-only fit cannot +separate domain size from beam coherence length; AFM, transverse scans, +rocking widths, or reciprocal-space maps are needed to constrain one of them. + +Constructing a ``PoissonSurface`` does **not** opt a calculation into +incoherent averaging. Wrapping the crystal is the explicit opt-in, and +existing scripts, saved crystals, and optimizer setups stay coherent:: + + from orgui.datautils.xrayutils.CTRincoherent import PoissonHeightDomains + from orgui.datautils.xrayutils.CTRopt import CTROptimizer + + # The crystal on its own is unchanged and fully coherent. + coherent = crystal.F2(h, k, l) + + domains = PoissonHeightDomains( + crystal, + surface="rough_surface", + incoherent_fraction=0.35, + ) + mixed = domains.F2(h, k, l) + +The fraction is a fixed setting until it is explicitly added as a fit +parameter, after which it occupies the first entry of the model block:: + + domains.addFitParameter( + "incoherent_fraction", + limits=(0.0, 1.0), + name="rough_surface incoherent_fraction", + ) + + fit = CTROptimizer(domains, ctrs) + fit.prepareFit() + assert fit.xtal is fit.model.coherent_model + +The wrapper is passed through the optimizer's existing model argument. +``optimizer.xtal`` remains the coherent crystal, so registered fit callbacks +and displacement constraints keep receiving an ``SXRDCrystal``, while +``optimizer.model`` is the fitted forward model. Resolution is applied to +``F2`` before the conversion back to a stored ``|F|``. + +The first implementation is kinematical only. Combining an incoherent +``F2`` model with DWBA raises during ``prepareFit``: changing the flat +surface height changes the optical reference profile and its internal fields, +so a correct DWBA ensemble would have to prepare and evaluate each height +separately and mix ``abs(r) ** 2``. The target surface must also be the +topmost component, since a component above it could have a position which +depends on the selected height. + +Height states are indexed by the structural layer :math:`n` of the top filled +layer, and the mass of that state is ``probability(n + 1)``: layer :math:`n` +is the top filled layer exactly when the signed height change equals +:math:`n + 1`. The retained interval is chosen from the calculated +probability masses and the profile's tail target, never from measured CTR +values, and is renormalized once. ``exact_layer_count`` (default 10) retains +every state for a narrow distribution; it is a policy switch, not a cap. + API reference ------------- @@ -702,9 +828,31 @@ API reference :member-order: bysource .. autoclass:: orgui.datautils.xrayutils.CTRcalc.SXRDCrystal - :members: F, F_surf, setGlobalReferenceUnitCell + :members: F, F2, F_surf, evaluate_kinematic, setGlobalReferenceUnitCell + :member-order: bysource + +.. autoclass:: orgui.datautils.xrayutils.CTRincoherent.IncoherentModel + :members: coherent_model, setParameters, validate, to_config + :member-order: bysource + +.. autoclass:: orgui.datautils.xrayutils.CTRincoherent.IncoherentF2Model + :members: F2 + :member-order: bysource + +.. autoclass:: orgui.datautils.xrayutils.CTRincoherent.CoherentStateEnsembleModel + :members: incoherent_fraction, accumulate_states :member-order: bysource +.. autoclass:: orgui.datautils.xrayutils.CTRincoherent.PoissonHeightDomains + :members: incoherent_fraction, target_surface + :member-order: bysource + +.. autofunction:: orgui.datautils.xrayutils.CTRincoherent.available_incoherent_models + +.. autofunction:: orgui.datautils.xrayutils.CTRincoherent.create_incoherent_model + +.. autofunction:: orgui.datautils.xrayutils.CTRincoherent.create_incoherent_model_from_config + .. autoclass:: orgui.datautils.xrayutils.CTRuc.UnitCell :members: F_uc, F_bulk, setReferenceUnitCell, supercell, affine_layer_transform, as_surface_termination :member-order: bysource @@ -718,7 +866,11 @@ API reference :member-order: bysource .. autoclass:: orgui.datautils.xrayutils.CTRfilm.PoissonSurface - :members: F_uc, uc_area, setReferenceUnitCell + :members: F_uc, uc_area, setReferenceUnitCell, flat_domain_corrections + :member-order: bysource + +.. autoclass:: orgui.datautils.xrayutils.CTRfilm.FlatHeightCorrections + :members: excluded_probability, iter_states, as_array :member-order: bysource .. autofunction:: orgui.datautils.xrayutils.CTRutil.generate_surface_termination_cells diff --git a/orgui/datautils/xrayutils/CTRincoherent.py b/orgui/datautils/xrayutils/CTRincoherent.py index 9e5a332..0b1ffef 100644 --- a/orgui/datautils/xrayutils/CTRincoherent.py +++ b/orgui/datautils/xrayutils/CTRincoherent.py @@ -423,8 +423,30 @@ def parameter_list(self): # -- composite serialization ---------------------------------------- + def _sync_local_parameter_values(self): + """Fill in local parameter values which were never set. + + ``addFitParameter`` records limits and a prior but leaves ``value`` + at ``None`` until the first ``setParameters`` call. Serializing that + state would store a parameter with no value, and restoring it raises + while rebuilding the basis. Reading the current values out of the + basis keeps the local subtree the single source for them. + """ + records = LinearFitFunctions.parameter_list(self) + if not records or all( + record.value is not None for record in records + ): + return + values = LinearFitFunctions.getStartParamAndLimits( + self, force_recalculate=True + )[0] + for record, value in zip(records, values): + if record.value is None: + record.value = float(value) + def _local_parameters_to_dict(self): """Return only the wrapper-local parameter subtree.""" + self._sync_local_parameter_values() return LinearFitFunctions.parametersToDict(self) def _local_parameters_from_dict(self, data, override_values=True): @@ -484,6 +506,15 @@ def to_config(self): them. The wrapped crystal is not included: the factory receives it separately, so the coherent structure keeps its own canonical file. + This dictionary is deliberately not part of any session or ``.xtal`` + format. A session layer may store it later and read an absent + dictionary as ordinary coherent behavior. + + What survives a round-trip is what the existing ``Parameter`` + serialization supports: names, limits, current values, errors, and + numeric priors. A prior held as a distribution object is dropped to + ``None`` by ``Parameter.asdict``; this feature does not widen that. + :rtype: dict """ if not self.model_type: @@ -831,10 +862,22 @@ def create_incoherent_model_from_config(crystal, config): record = _registered(config["type"]) settings = config.get("settings") or {} - model = record.model_class(crystal, **dict(settings)) + if not isinstance(settings, dict): + raise ValueError("Incoherent model settings must be a dict") + try: + model = record.model_class(crystal, **dict(settings)) + except TypeError as error: + raise ValueError( + f"Settings {sorted(settings)} do not construct " + f"{config['type']!r}: {error}" + ) from error + parameters = config.get("parameters") - if parameters: - model._local_parameters_from_dict(parameters) + if parameters is None: + return model + if not isinstance(parameters, dict): + raise ValueError("Incoherent model parameters must be a dict") + model._local_parameters_from_dict(parameters) return model diff --git a/orgui/datautils/xrayutils/test/test_CTRincoherent.py b/orgui/datautils/xrayutils/test/test_CTRincoherent.py index ad1f9e3..9b48a38 100644 --- a/orgui/datautils/xrayutils/test/test_CTRincoherent.py +++ b/orgui/datautils/xrayutils/test/test_CTRincoherent.py @@ -1021,5 +1021,116 @@ def test_survives_the_deep_copy_the_optimizer_performs(self): ) +class TestConfigurationRoundTrip(PoissonModelMixin): + """What a dictionary round-trip preserves, and what it rejects.""" + + PROFILE = PoissonProfile(1.5, alpha=0.5, tail_probability=1e-14) + + def fitted_model(self): + """Return a model whose fraction is fitted, bounded, and errored.""" + model, crystal, _, _ = self.model(self.PROFILE, kappa=0.35) + model.addFitParameter( + "incoherent_fraction", + limits=(0.1, 0.9), + name="surface incoherent_fraction", + prior=(0.2, 0.8), + ) + model.setFitErrors( + np.concatenate( + ([0.05], np.zeros(model.n_coherent_parameters)) + ) + ) + return model, crystal + + def test_round_trip_preserves_the_local_block(self): + """Type, settings, value, limits, errors, and priors all survive.""" + model, crystal = self.fitted_model() + config = model.to_config() + + restored = create_incoherent_model_from_config(crystal, config) + self.assertEqual(restored.model_type, model.model_type) + self.assertEqual(restored.surface, model.surface) + self.assertEqual( + restored.exact_layer_count, model.exact_layer_count + ) + self.assertEqual(restored.fitparnames, model.fitparnames) + np.testing.assert_allclose( + restored.getInitialParameters(), model.getInitialParameters() + ) + + _, lower, upper = restored.getStartParamAndLimits() + self.assertAlmostEqual(lower[0], 0.1) + self.assertAlmostEqual(upper[0], 0.9) + np.testing.assert_allclose(restored.priors[0], (0.2, 0.8)) + np.testing.assert_allclose(restored.getFitErrors()[0], 0.05) + np.testing.assert_allclose( + restored.F2(self.H, self.K, self.L), + model.F2(self.H, self.K, self.L), + rtol=1e-14, + ) + + def test_the_fraction_is_stored_once(self): + """Settings must not duplicate the local basis value.""" + model, _ = self.fitted_model() + config = model.to_config() + self.assertNotIn("incoherent_fraction", config["settings"]) + np.testing.assert_allclose( + config["parameters"]["basis_0"], [0.35] + ) + + def test_a_fixed_fraction_round_trips_without_a_fit_parameter(self): + """A setting stays fixed and still survives the round-trip.""" + model, crystal, _, _ = self.model(self.PROFILE, kappa=0.6) + self.assertEqual(model.n_local_parameters, 0) + restored = create_incoherent_model_from_config( + crystal, model.to_config() + ) + self.assertAlmostEqual(restored.incoherent_fraction, 0.6) + self.assertEqual(restored.n_local_parameters, 0) + + def test_missing_settings_fail_with_the_offending_keys(self): + """A configuration which cannot construct its model says so.""" + model, crystal, _, _ = self.model(self.PROFILE) + config = model.to_config() + del config["settings"]["surface"] + with self.assertRaisesRegex(ValueError, "do not construct"): + create_incoherent_model_from_config(crystal, config) + + config = model.to_config() + config["settings"]["nonsense"] = 1 + with self.assertRaisesRegex(ValueError, "do not construct"): + create_incoherent_model_from_config(crystal, config) + + def test_malformed_configurations_are_rejected(self): + """Deserialization never guesses at a broken dictionary.""" + model, crystal, _, _ = self.model(self.PROFILE) + for bad, pattern in ( + ("not a dict", "must be a dict"), + ({}, "needs a type"), + ({"type": "poisson_height_domains", "settings": 5}, "must be a dict"), + ): + with self.subTest(config=bad): + with self.assertRaisesRegex(ValueError, pattern): + create_incoherent_model_from_config(crystal, bad) + + config = model.to_config() + config["parameters"] = ["not", "a", "dict"] + with self.assertRaisesRegex(ValueError, "must be a dict"): + create_incoherent_model_from_config(crystal, config) + + def test_the_wrapped_crystal_is_supplied_not_stored(self): + """The coherent structure keeps its own canonical file.""" + model, _, _, _ = self.model(self.PROFILE, kappa=0.4) + config = model.to_config() + self.assertEqual(set(config), {"type", "settings", "parameters"}) + self.assertNotIn("coherent", config["parameters"]) + + _, other_crystal, _, _ = self.model(self.PROFILE) + restored = create_incoherent_model_from_config( + other_crystal, config + ) + self.assertIs(restored.coherent_model, other_crystal) + + if __name__ == "__main__": unittest.main() From 06b4e074e082db600e05961243a790756c640cb1 Mon Sep 17 00:00:00 2001 From: "T. Fuchs" Date: Sat, 12 Sep 2026 14:38:10 -0400 Subject: [PATCH 13/13] feat(phys): allow components stacked above a rough surface An overlayer is placed once at the mean surface height and never moves while the height states stream, so it belongs in the common amplitude exactly as the coherent model already treats it. --- CHANGELOG.md | 6 +- doc/design/incoherent_ctr_models.md | 41 ++++++---- doc/source/ctr_structure_factors.rst | 12 ++- orgui/datautils/xrayutils/CTRincoherent.py | 45 ++++++----- .../xrayutils/test/test_CTRincoherent.py | 76 +++++++++++++++++-- 5 files changed, 133 insertions(+), 47 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 506db3a..4b87c62 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,8 +25,10 @@ Scientific and analysis additions: `SXRDCrystal` gains `F2` and `evaluate_kinematic`, resolution now acts on `F2` before the conversion back to a stored amplitude, and `CTROptimizer` gains `n_parameters` for the full prepared vector. The first implementation - is kinematical: combining it with DWBA, or targeting a surface which is not - the topmost component, fails explicitly during `prepareFit`. + is kinematical: combining it with DWBA fails explicitly during + `prepareFit`. A component such as a water layer may be stacked above the + target surface, where it is common to every domain at the mean surface + height, as in the coherent model. - **Added live DWBA predictions to CTR fitting.** ``CTROptimizer.set_dwba`` now evaluates the optimizer-owned crystal through the semi-infinite DWBA diff --git a/doc/design/incoherent_ctr_models.md b/doc/design/incoherent_ctr_models.md index e945bd2..cedacd1 100644 --- a/doc/design/incoherent_ctr_models.md +++ b/doc/design/incoherent_ctr_models.md @@ -1003,14 +1003,28 @@ Object identity or only `_basis_created` is not an adequate cache key. ### Stacking boundary -The first implementation targets one `PoissonSurface` which is the highest -material component in stacking order. This is the case in which all material -below the surface is common and only the surface correction varies. - -If a component is stacked above the target surface, its absolute position or -structure may depend on the selected height. Treating it as part of -`A_common` would be wrong. Initially reject that topology with an actionable -error. A later extension can split the evaluation into: +The target `PoissonSurface` does not have to be the highest material +component. This record originally required it to be, reasoning that a +component stacked above it "may" depend on the selected height and so could +not join `A_common`. Checking the implementation showed that dependence does +not exist: `apply_stacking` places an overlayer at the surface's +`stacking_height_absolute`, which is `mean_height_absolute`, a function of +the fit parameters `W` and `offset` and not of any height state. Nothing +re-stacks while the states are streamed, so an overlayer holds one position +for every domain. Measured directly, a Film stacked on a rough surface keeps +one `below_H` and one `F_uc` across all streamed states. + +An overlayer therefore belongs in `A_common`, and putting it there reproduces +exactly what the coherent model already does with it. The `kappa = 0` +endpoint stays bitwise equal to `SXRDCrystal.F2`. + +What that costs is a physical approximation worth stating: the overlayer sits +at the *mean* surface height in every domain rather than following each +domain's own height. In a strict large-domain limit, water or a cap would +ride each terrace. Reproducing that is the deferred extension, and it is a +change of physics rather than a lifted restriction, because it would make the +incoherent model more correct than the coherent baseline and so break the +`kappa = 0` identity. It would split the evaluation into: - a cached common prefix below the target surface; - the target flat-height correction; @@ -1376,11 +1390,10 @@ Work: Add the positive-integer non-fit setting `exact_layer_count`, default 10. 2. Accept the coherent crystal as the first positional argument and a stable target surface name. Validate exactly one matching `PoissonSurface`, a - finite fraction in `[0, 1]`, and the supported stacking topology. For the - first implementation, the target must be the final component in resolved - stacking order, must already be bound immediately above its Film, and may - have no same-level or later component whose position could depend on the - selected height. + finite fraction in `[0, 1]`, and that the target is bound immediately + above its Film. The target need not be the final component in resolved + stacking order: see "Stacking boundary" for why an overlayer is common to + every domain rather than state-dependent. 3. Use the context to evaluate the common crystal amplitude once and to turn raw flat-height corrections into complete state amplitudes with all area, weight, and outer-domain factors applied. @@ -1614,7 +1627,7 @@ The first implementation is complete only when: - the optimizer consumes the wrapper through its existing model/parameter contract; - resolution acts on `F2` before conversion to `|F|`; -- unsupported DWBA and stacking cases fail explicitly; +- unsupported DWBA and non-Poisson targets fail explicitly; - bulk and common Film amplitudes are evaluated once per HKL batch, verified by call-count instrumentation as well as timing; and - no long-lived cache with incomplete invalidation is introduced. diff --git a/doc/source/ctr_structure_factors.rst b/doc/source/ctr_structure_factors.rst index 2c38a84..21a7960 100644 --- a/doc/source/ctr_structure_factors.rst +++ b/doc/source/ctr_structure_factors.rst @@ -798,9 +798,15 @@ The first implementation is kinematical only. Combining an incoherent ``F2`` model with DWBA raises during ``prepareFit``: changing the flat surface height changes the optical reference profile and its internal fields, so a correct DWBA ensemble would have to prepare and evaluate each height -separately and mix ``abs(r) ** 2``. The target surface must also be the -topmost component, since a component above it could have a position which -depends on the selected height. +separately and mix ``abs(r) ** 2``. + +The target surface does not have to be the topmost component: a water layer +or a cap may be stacked above it. Anything above the surface is placed once +at the surface's *mean* height and is common to every domain, which is how +the coherent model already treats it, so the :math:`\kappa = 0` endpoint is +unchanged. Note the approximation this carries: in a strict large-domain +limit an overlayer would follow each domain's own height rather than the +mean. Height states are indexed by the structural layer :math:`n` of the top filled layer, and the mass of that state is ``probability(n + 1)``: layer :math:`n` diff --git a/orgui/datautils/xrayutils/CTRincoherent.py b/orgui/datautils/xrayutils/CTRincoherent.py index 0b1ffef..c19c6fc 100644 --- a/orgui/datautils/xrayutils/CTRincoherent.py +++ b/orgui/datautils/xrayutils/CTRincoherent.py @@ -907,8 +907,8 @@ class PoissonHeightDomains(CoherentStateEnsembleModel): :param str name: Wrapper name, the default prefix for its local parameter. :raises ValueError: - If the target is missing, duplicated, not a ``PoissonSurface``, not - the topmost component, or the fraction is out of range. + If the target is missing, duplicated, not a ``PoissonSurface``, or + the fraction is out of range. """ model_type = "poisson_height_domains" @@ -992,7 +992,17 @@ def target_surface(self): return matches[0] def _validate_target(self): - """Check the target type and the supported stacking topology.""" + """Check that the named component is a Poisson surface. + + The target does not have to be the topmost component. Anything + stacked above it is placed by ``apply_stacking`` at the surface's + ``stacking_height_absolute``, which is the *mean* surface height: a + function of the fit parameters, not of the height state being + evaluated. Nothing re-stacks while the states are streamed, so an + overlayer holds one position for every domain and belongs in the + common amplitude, exactly as the coherent model already treats it. + See :meth:`_iter_states` for what that approximation costs. + """ target = self.target_surface() if not isinstance(target, PoissonSurface): raise ValueError( @@ -1000,23 +1010,6 @@ def _validate_target(self): f"{type(target).__name__}, not a PoissonSurface" ) - crystal = self._coherent_model - ordered = list(crystal.uc_surface_list_ordered) - if not ordered or ordered[-1] is not target: - raise ValueError( - f"{self.surface!r} must be the topmost component in stacking " - "order. A component above it could have a position which " - "depends on the selected height, so it cannot be treated as " - "part of the common amplitude" - ) - levels = np.asarray(crystal.uc_stacking_ordered) - if levels.size and int(np.sum(levels == levels[-1])) != 1: - raise ValueError( - f"{self.surface!r} shares its stacking level with another " - "component, whose position could depend on the selected " - "height" - ) - def _validate_model(self, forward_model): """Check that the target is bound to the Film it corrects.""" target = self._bound_target() @@ -1034,7 +1027,17 @@ def _bound_target(self): return self.target_surface() def _iter_states(self, context): - """Stream complete crystal amplitudes, one flat height at a time.""" + """Stream complete crystal amplitudes, one flat height at a time. + + Every component other than the target is common to all states, + including any stacked above the surface. Such an overlayer therefore + sits at the mean surface height in every domain rather than following + each domain's own height. That is the same approximation the coherent + model makes, so the ``kappa = 0`` endpoint still reproduces + ``SXRDCrystal.F2`` exactly; it is not the large-domain limit an + overlayer would strictly have, which would require re-stacking and + re-evaluating the suffix for each height. + """ def state_evaluator(component, h, k, l): # noqa: E741 return component.flat_domain_corrections( diff --git a/orgui/datautils/xrayutils/test/test_CTRincoherent.py b/orgui/datautils/xrayutils/test/test_CTRincoherent.py index 9b48a38..3554807 100644 --- a/orgui/datautils/xrayutils/test/test_CTRincoherent.py +++ b/orgui/datautils/xrayutils/test/test_CTRincoherent.py @@ -13,7 +13,7 @@ import numpy as np -from .. import CTRcalc, CTRincoherent +from .. import CTRcalc, CTRfilm, CTRincoherent from ..CTRincoherent import ( CoherentStateEnsembleModel, IncoherentF2Model, @@ -971,14 +971,14 @@ def test_rejects_a_target_which_is_not_a_poisson_surface(self): with self.assertRaisesRegex(ValueError, "not a PoissonSurface"): CTRincoherent.PoissonHeightDomains(crystal, surface="film") - def test_rejects_a_component_stacked_above_the_target(self): - """A component above the surface may move with the chosen height.""" + def capped_crystal(self): + """Return a crystal with an overlayer stacked above the surface.""" w_base = _poisson_oracle.minimum_base_width(self.PROFILE) crystal, surface = _poisson_oracle.poisson_crystal( self.PROFILE, w_base=w_base ) - cap = CTRcalc.UnitCell([3.0, 3.0, 6.0], [90.0, 90.0, 90.0], name="cap") - cap.addAtom("O", [0.0, 0.0, 0.0], 0.1, 0.1, 1.0, layer=0) + cap = CTRfilm.Film(_poisson_oracle.layered_cell(2, "cap"), name="cap") + cap.basis[0] = 2.0 stacked = CTRcalc.SXRDCrystal( _poisson_oracle.layered_cell(2, "bulk"), crystal.uc_surface_list[0], @@ -986,8 +986,70 @@ def test_rejects_a_component_stacked_above_the_target(self): cap, stacking=np.array([1, 2, 3]), ) - with self.assertRaisesRegex(ValueError, "topmost component"): - CTRincoherent.PoissonHeightDomains(stacked, surface="surface") + stacked.apply_stacking() + return stacked, surface, cap + + def test_accepts_a_component_stacked_above_the_target(self): + """An overlayer does not have to be refused. + + Anything above the surface is placed once by ``apply_stacking`` at + the surface's mean height, and nothing re-stacks while the states are + streamed. Such a component therefore holds one position for every + domain and belongs in the common amplitude, which is exactly how the + coherent model already treats it. + """ + stacked, _, _ = self.capped_crystal() + self.assertIsNot(stacked.uc_surface_list_ordered[-1].name, "surface") + + model = CTRincoherent.PoissonHeightDomains( + stacked, surface="surface", incoherent_fraction=0.0 + ) + # The coherent endpoint still reproduces the crystal exactly. + np.testing.assert_array_equal( + model.F2(self.H, self.K, self.L), + stacked.F2(self.H, self.K, self.L), + ) + + mixed = CTRincoherent.PoissonHeightDomains( + stacked, surface="surface", incoherent_fraction=1.0 + ) + self.assertGreater( + np.max( + np.abs( + mixed.F2(self.H, self.K, self.L) + - model.F2(self.H, self.K, self.L) + ) + ), + 0.0, + ) + + def test_an_overlayer_is_common_to_every_state(self): + """The overlayer contributes identically to each height state.""" + stacked, _, cap = self.capped_crystal() + context = CTRincoherent.KinematicIncoherentContext( + stacked, self.H, self.K, self.L + ) + cap_amplitude = context.component("cap").amplitude + + # It sits inside the common amplitude, not inside any state. + common = context.common_amplitude("surface") + np.testing.assert_allclose( + common, + context.coherent.bulk + + context.component("film").amplitude + + cap_amplitude, + rtol=1e-12, + ) + + # And its own placement does not move while the states stream. + model = CTRincoherent.PoissonHeightDomains( + stacked, surface="surface", incoherent_fraction=1.0 + ) + positions = { + round(float(cap.below_H), 9) + for _ in model._iter_states(context) + } + self.assertEqual(len(positions), 1) def test_rejects_an_invalid_state_count(self): """``exact_layer_count`` is a positive integer policy setting."""