diff --git a/CHANGELOG.md b/CHANGELOG.md index d239056..4b87c62 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,29 @@ 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 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 model, forms predictions independently in each dataset's stored F or diff --git a/doc/design/incoherent_ctr_models.md b/doc/design/incoherent_ctr_models.md new file mode 100644 index 0000000..cedacd1 --- /dev/null +++ b/doc/design/incoherent_ctr_models.md @@ -0,0 +1,1711 @@ +# Incoherent CTR models: design and implementation plan + +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 +`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. + +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, `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 +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 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 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 +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 setFitParameters(self, values): # raises; use setParameters + ... + 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", +) +``` + +`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 +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. + +`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. + +`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", +) +# 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) +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 model-block snapshot; the full + # vector is get_parameters() +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 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. + +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 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 +nonzero states, retain them all. For a wider support, retain a contiguous +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. + +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 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; +- 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`. 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 + 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. + +### 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` +- 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: + +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 + +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 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 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 + 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. 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`. +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. 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 + 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. +- 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. +- 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. 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. + 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. 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 + 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 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 + 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. + 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, 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. +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, 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. +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`, `bounds`, `fitparnames`, `get_parameters()` and + `n_parameters` describe the same full vector with every optional prefix + 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 + 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. + +### 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 +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_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 +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 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. + +## 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. diff --git a/doc/source/ctr_structure_factors.rst b/doc/source/ctr_structure_factors.rst index dbb9c47..21a7960 100644 --- a/doc/source/ctr_structure_factors.rst +++ b/doc/source/ctr_structure_factors.rst @@ -684,6 +684,138 @@ 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 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` +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 +834,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 +872,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/CTRcalc.py b/orgui/datautils/xrayutils/CTRcalc.py index 38ab076..25ae1a3 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,133 @@ 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] + 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]: + yield matrix, area_scale * occup * weight + + 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/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/CTRincoherent.py b/orgui/datautils/xrayutils/CTRincoherent.py new file mode 100644 index 0000000..c19c6fc --- /dev/null +++ b/orgui/datautils/xrayutils/CTRincoherent.py @@ -0,0 +1,1054 @@ +# /*########################################################################## +# +# 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 .CTRfilm import PoissonSurface +from .CTRutil import LinearFitFunctions + +__all__ = [ + "CoherentStateEnsembleModel", + "IncoherentF2Model", + "IncoherentModel", + "IncoherentModelInfo", + "IncoherentState", + "KinematicIncoherentContext", + "PoissonHeightDomains", + "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 _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): + """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. + + 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: + 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 {} + 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 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 + + +@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``, 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 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( + f"Component {self.surface!r} is a " + f"{type(target).__name__}, not a PoissonSurface" + ) + + 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. + + 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( + 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/CTRopt.py b/orgui/datautils/xrayutils/CTRopt.py index fa1c781..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 @@ -1412,11 +1518,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/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/_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..97ababf 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,439 @@ 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 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"): 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() diff --git a/orgui/datautils/xrayutils/test/test_CTRincoherent.py b/orgui/datautils/xrayutils/test/test_CTRincoherent.py new file mode 100644 index 0000000..3554807 --- /dev/null +++ b/orgui/datautils/xrayutils/test/test_CTRincoherent.py @@ -0,0 +1,1198 @@ +"""Tests for the incoherent CTR wrapper and the Poisson height-domain model. + +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 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 +import unittest + +import numpy as np + +from .. import CTRcalc, CTRfilm, 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 ..CTRdistributions import PoissonProfile +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)) + + +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 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 = 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], + surface, + cap, + stacking=np.array([1, 2, 3]), + ) + 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.""" + 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, + ) + + +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() diff --git a/orgui/datautils/xrayutils/test/test_CTRopt.py b/orgui/datautils/xrayutils/test/test_CTRopt.py index 06d7f1c..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): @@ -1636,6 +1645,468 @@ 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 _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()