diff --git a/docs/plans/nrss_q_coordinate_semantics.md b/docs/plans/nrss_q_coordinate_semantics.md new file mode 100644 index 00000000..9564016a --- /dev/null +++ b/docs/plans/nrss_q_coordinate_semantics.md @@ -0,0 +1,573 @@ +# NRSS q-coordinate semantics unification plan + +## Purpose + +This note records the implementation and validation plan for unifying q-coordinate semantics between: + +- NRSS simulation outputs +- PyHyperScattering reduction workflows +- downstream analytical form-factor comparisons + +It is intentionally stored outside the Sphinx tree so it can be used as a dev-facing design and execution note without changing docs tooling. + +## Current diagnosis + +The mismatch described in `/homes/deand/dev/NRSS/test-reports/nrss_q_coordinate_reduction_issue_draft.md` is real and is visible in current code: + +- `WPIntegrator` in PyHyperScattering is a warp-polar remesher over existing `qx/qy` coordinates. + - It finds the image center from `qx=0` and `qy=0`. + - It runs `warp_polar`. + - It labels the radial axis with `sqrt(qx^2 + qy^2)`. + - This is `q_perp`, not a detector-geometry-corrected `|q|`. +- NRSS and `cyrsoxsLoader` currently export FFT-style `qx/qy` detector-plane coordinates derived from `PhysSize` and array shape. +- The modern NRSS 3D backend already uses detector-aware projection geometry internally during intensity formation, but that geometry is not preserved in the exported xarray metadata. +- As a result: + - 2D NRSS reduced with `WPIntegrator` is semantically fine. + - 3D NRSS reduced with `WPIntegrator` carries a radial axis that is semantically wrong for comparison to geometry-aware experiment reduction or pure analytical `I(|q|)`. + +## Repository state at time of writing + +### PyHyperScattering + +- Repo: `/homes/deand/dev/PyHyperScattering` +- Active branch: + - `221-unify-q-coordinate-semantics-between-nrss-simulation-outputs-and-pyhyperscattering-reduction-workflows` +- Worktree state: + - clean +- Important finding: + - this branch is currently identical to `main` + - there is no implementation on the branch yet + +### NRSS + +- Repo: `/homes/deand/dev/NRSS` +- Active branch: + - `backend/cupy-cyrsoxs-clone` +- Worktree state: + - dirty + - modified: `src/NRSS/backends/cupy_rsoxs.py` + - several deleted design/spec markdown files in repo root +- Consequence: + - phase 1 should avoid requiring NRSS code changes + - local NRSS should initially be used only as a consumer/validation target + +## Key code anchors + +### PyHyperScattering + +- `WPIntegrator` + - `/homes/deand/dev/PyHyperScattering/src/PyHyperScattering/WPIntegrator.py` +- `cyrsoxsLoader` + - `/homes/deand/dev/PyHyperScattering/src/PyHyperScattering/cyrsoxsLoader.py` +- `PFGeneralIntegrator` + - `/homes/deand/dev/PyHyperScattering/src/PyHyperScattering/PFGeneralIntegrator.py` +- Existing tests + - `/homes/deand/dev/PyHyperScattering/tests/test_WPintegrator.py` + - `/homes/deand/dev/PyHyperScattering/tests/test_pf_general_integrator_basic.py` + +### NRSS + +- legacy xarray export + - `/homes/deand/dev/NRSS/src/NRSS/morphology.py` +- cupy backend xarray export + - `/homes/deand/dev/NRSS/src/NRSS/backends/cupy_rsoxs.py` +- maintained sphere form-factor validation + - `/homes/deand/dev/NRSS/tests/validation/test_analytical_sphere_form_factor.py` +- existing plot-writing validation pattern + - `/homes/deand/dev/NRSS/tests/validation/test_sphere_contrast_scaling.py` + - `/homes/deand/dev/NRSS/tests/validation/test_sphere_orientational_contrast_scaling.py` + +## High-level decision + +Do **not** force the first implementation through `PFGeneralIntegrator` or synthetic pyFAI metadata. + +Reason: + +- `PFGeneralIntegrator` is built around raw detector images with `pix_x/pix_y` plus pyFAI-style geometry metadata. +- NRSS currently provides `energy/qy/qx` arrays, not raw detector pixel stacks with a stable calibration contract. +- inventing fake pyFAI geometry in phase 1 adds risk and does not solve the metadata problem + +Therefore: + +- build a new PyHyperScattering-side `NRSSIntegrator` first +- keep it native to NRSS-style `qx/qy` arrays +- use explicit semantics and attrs +- validate with local NRSS before any NRSS metadata bridge is merged + +## Planned end state + +### For 2D NRSS output + +- reduction remains a radial remesh over detector-plane reciprocal coordinates +- the semantics should be explicit: + - radial axis is `q_perp` +- the result may still use a dimension named `q` for compatibility, but must carry metadata that says it is `q_perp` + +### For 3D NRSS output + +- reduction should still use the simulated detector-plane image intensity +- but the radial coordinate should be corrected from detector-plane `q_perp` to full `|q|` +- this should match the geometry already used inside the 3D NRSS detector-projection backend + +### Migration target + +- short term: + - keep detector-space `WPIntegrator`-based validations where they already exist + - add `NRSSIntegrator` validations against pure analytical form factors +- long term: + - replace detector-space analytical comparison tests in NRSS with `NRSSIntegrator` vs pure analytical `I(|q|)` tests + +## Flow 1: PyHyperScattering implementation first + +This is the required first step. + +### Scope + +Implement a new `NRSSIntegrator` class in PyHyperScattering. + +Suggested file: + +- `/homes/deand/dev/PyHyperScattering/src/PyHyperScattering/NRSSIntegrator.py` + +Expected exports also need to be wired through whatever import surface currently exposes integrators, for example: + +- `src/PyHyperScattering/integrate.py` +- package `__init__` if needed by current conventions + +### Required behavior + +#### Input contract + +`NRSSIntegrator` should accept an xarray with NRSS-like coordinates: + +- spatial dims: + - `qy`, `qx` +- optional stack dims: + - `energy` + - any additional grouping/index coords that current integrators already tolerate + +It should also accept explicit metadata kwargs for cases where attrs are missing: + +- `phys_size_nm` +- `shape_zyx` or at least `z_dim` +- `energy_ev` or an energy coord +- `projection_mode` or equivalent semantic hint + +It should first look in `img.attrs`, then use explicit kwargs as fallback. + +### Semantic mode detection + +The integrator must distinguish at least two modes: + +#### Mode A: reciprocal-plane / 2D semantics + +Use this when: + +- `z_dim == 1` +- or attrs explicitly say the output is reciprocal-plane / 2D + +Behavior: + +- reduce exactly like `WPIntegrator` +- radial axis semantics = `q_perp` + +#### Mode B: detector-aware / 3D semantics + +Use this when: + +- `z_dim > 1` +- or attrs explicitly say detector-aware / 3D projection + +Behavior: + +- use the same detector-plane intensity image +- keep the same angular remesh approach +- but relabel the radial axis using corrected `|q|` + +### Core geometry formula for 3D mode + +This formula should match the logic already present in NRSS cupy backend detector-projection geometry. + +Definitions: + +- `q_perp = sqrt(qx^2 + qy^2)` +- `k = 2*pi / lambda` +- with energy in eV: + - `lambda_nm = 1239.84197 / energy_ev` + - `k = 2*pi / lambda_nm` + +Detector-aware projection z component: + +```text +qz = -k + sqrt(k^2 - q_perp^2) +``` + +Corrected magnitude: + +```text +q = sqrt(q_perp^2 + qz^2) +``` + +Constraints: + +- values are only valid where `k^2 - q_perp^2 >= 0` +- if the detector-plane image already contains invalid edge regions as NaN, preserve that masking behavior + +### Output contract + +Return a `chi/q` reduced xarray compatible with existing downstream usage. + +Recommended output shape: + +- `chi`, `q` +- plus any stacked axes preserved the same way `WPIntegrator` preserves them + +Required metadata: + +- `radial_semantics` + - `"q_perp"` for 2D mode + - `"q_abs_detector_corrected"` for 3D mode +- `source_integrator` + - `"NRSSIntegrator"` +- `nrss_semantic_mode` + - `"2d_reciprocal_plane"` or `"3d_detector_aware"` +- `phys_size_nm` +- `z_dim` +- `energy_ev` if scalar + +Recommended auxiliary coordinate for 3D mode: + +- keep output dim name `q` for compatibility +- also attach a `q_perp` coordinate or attr summary when possible + +### Compatibility rule + +Do **not** change `WPIntegrator` behavior in phase 1. + +Only add a clarifying docstring or warning so users understand: + +- `WPIntegrator` radial output is `q_perp` +- it is still correct for 2D NRSS reciprocal-plane output +- it is not the geometry-correct reduction for 3D detector-aware NRSS output + +### PyHyperScattering tests to add + +Add a new test module: + +- `/homes/deand/dev/PyHyperScattering/tests/test_NRSSIntegrator.py` + +At minimum include: + +1. `test_nrss_integrator_2d_matches_wp_semantics` + - build a synthetic `qx/qy` image + - set `z_dim=1` + - verify reduced radial axis equals `WPIntegrator` radial axis + - verify attrs say `radial_semantics == "q_perp"` + +2. `test_nrss_integrator_3d_corrects_radial_axis` + - build a synthetic `qx/qy` image with an energy + - set `z_dim > 1` + - verify output `q` differs from `sqrt(qx^2 + qy^2)` at higher q + - verify corrected values match the explicit formula above + +3. `test_nrss_integrator_falls_back_to_explicit_kwargs` + - no attrs on the xarray + - pass explicit `phys_size_nm`, `z_dim`, and `energy_ev` + - verify reduction succeeds + +4. `test_nrss_integrator_preserves_stack_axis` + - use an `energy` stack or another grouping coord + - verify output preserves group axis + +5. `test_wp_integrator_semantics_note` + - if behavior is not changed, at least assert the warning/doc metadata path if one is added + +### PyHyperScattering docs to update + +Because this plan is outside Sphinx, no docs build change is required now. + +But after implementation, update at least: + +- `/homes/deand/dev/PyHyperScattering/docs/source/getting_started/integration.rst` + +Add short text clarifying: + +- `WPIntegrator` is a `qx/qy -> chi/q_perp` polar remesher +- `NRSSIntegrator` is the correct reducer for NRSS detector-aware 3D outputs + +## Flow 2: NRSS validation and migration + +This phase starts after `NRSSIntegrator` exists and has passing unit tests in PyHyperScattering. + +### Immediate goal + +Use local NRSS as the validation target for `NRSSIntegrator` without requiring NRSS metadata changes first. + +### Environment notes + +Per local instructions: + +- for joint NRSS + PyHyperScattering work use environment: + - `nrss-pyhyper-dev` +- prefer `mamba` for environment work +- do not fire more than one env command at once + +### First validation target + +Base the dev-only validation style on: + +- `/homes/deand/dev/NRSS/tests/validation/test_analytical_sphere_form_factor.py` + +This file already provides: + +- sphere morphology construction +- analytical sphere form-factor comparison +- backend runtime path helpers +- optional plot writing to `test-reports/` + +Existing test-report plot pattern is also visible in: + +- `/homes/deand/dev/NRSS/tests/validation/test_sphere_contrast_scaling.py` +- `/homes/deand/dev/NRSS/tests/validation/test_sphere_orientational_contrast_scaling.py` + +They use: + +- `PLOT_DIR = REPO_ROOT / "test-reports" / "-dev"` +- `WRITE_VALIDATION_PLOTS = os.environ.get("NRSS_WRITE_VALIDATION_PLOTS", "").strip() == "1"` + +Follow that exact pattern. + +### NRSS dev-only tests to add + +These should initially be dev validation tests, not immediately promoted to maintained reference parity until the semantics settle. + +Suggested new test module: + +- `/homes/deand/dev/NRSS/tests/validation/test_nrss_integrator_sphere_form_factor.py` + +If the team prefers to keep it explicitly provisional, use a clearer name such as: + +- `/homes/deand/dev/NRSS/tests/validation/test_nrss_integrator_sphere_form_factor_dev.py` + +The test should: + +1. build or reuse the same sphere morphology style as `test_analytical_sphere_form_factor.py` +2. run NRSS backend(s) to produce detector-plane xarray output +3. reduce that output twice: + - once with `WPIntegrator` + - once with `NRSSIntegrator` +4. compare both reductions to the same pure analytical sphere form factor +5. write plots when `NRSS_WRITE_VALIDATION_PLOTS=1` + +### Required visualizations + +Each validation run should generate a comparison figure showing at minimum: + +1. detector image in `qx/qy` +2. `WPIntegrator` reduced curve vs detector-space expectation if still useful +3. `NRSSIntegrator` reduced curve vs pure analytical `I(|q|)` +4. residual or log-error plot +5. minima alignment markers if available + +The important expected result is: + +- `NRSSIntegrator` should align with the pure analytical form factor +- `WPIntegrator` should visibly diverge in the 3D case at higher q + +### Suggested plot/report directory + +Use a new dev report directory such as: + +- `/homes/deand/dev/NRSS/test-reports/nrss-integrator-sphere-dev` + +### Suggested scenario matrix + +At minimum: + +#### 2D scenario + +- `z_dim == 1` +- expected: + - `WPIntegrator` and `NRSSIntegrator` should agree + - both should align with the 2D-appropriate analytical interpretation + +#### 3D scenario + +- `z_dim > 1` +- expected: + - `WPIntegrator` follows detector-plane `q_perp` + - `NRSSIntegrator` follows analytical `|q|` + - `NRSSIntegrator` should provide the physically correct comparison target + +Possible concrete settings: + +- reuse sphere diameters already maintained: + - `70 nm` + - `128 nm` +- start with one energy: + - `285.0 eV` +- reuse `PhysSize` and shape conventions already established in the analytical sphere test where practical + +### Threshold strategy + +For the first NRSS dev-only validation, do not start with aggressive fixed thresholds. + +Instead: + +1. generate plots +2. measure: + - RMS log error + - P95 absolute log error + - minima position error + - integrated intensity agreement over the chosen q window +3. inspect results +4. then set thresholds appropriate for maintained migration + +### Long-term NRSS migration + +After `NRSSIntegrator` is validated: + +1. update NRSS maintained analytical sphere tests so the primary analytical comparison path is: + - `NRSS output -> NRSSIntegrator -> pure analytical I(|q|)` +2. demote or remove detector-space `WPIntegrator` analytical checks where they are only preserving the old semantics +3. keep a small compatibility test proving that: + - `WPIntegrator` still behaves as detector-plane `q_perp` + - it remains correct for 2D reciprocal-plane reductions + +## Metadata bridge to add later in NRSS + +This is intentionally **not** part of flow 1. + +After PyHyperScattering-side reduction is working and validated, NRSS should export enough attrs for automatic mode detection. + +### Suggested attrs on exported xarray + +- `phys_size_nm` +- `num_zyx` +- `z_dim` +- `nrss_output_semantics` + - `"2d_reciprocal_plane"` + - `"3d_detector_aware"` +- `backend` +- `backend_options` +- `energy_ev` if scalar or rely on energy coord + +### Likely files to patch later + +- `/homes/deand/dev/NRSS/src/NRSS/morphology.py` +- `/homes/deand/dev/NRSS/src/NRSS/backends/cupy_rsoxs.py` + +### Important caution + +Do not mix the metadata-bridge change into the first PyHyperScattering implementation unless absolutely necessary. + +Reason: + +- PyHyper should first be able to operate on current local NRSS output via explicit kwargs +- that keeps the semantic correction testable before any cross-repo contract change lands + +## Concrete execution order + +### Step 1 + +In PyHyperScattering: + +- add `NRSSIntegrator` +- add unit tests +- keep `WPIntegrator` behavior unchanged + +### Step 2 + +Smoke-test `NRSSIntegrator` on synthetic xarrays in PyHyperScattering only. + +### Step 3 + +Use local NRSS output directly with `NRSSIntegrator` by passing explicit metadata kwargs. + +### Step 4 + +Add NRSS dev-only sphere-form-factor validations with optional plot writing. + +### Step 5 + +Confirm expected semantic split: + +- 2D: + - `WPIntegrator == NRSSIntegrator` +- 3D: + - `NRSSIntegrator` matches pure analytical form factor better than `WPIntegrator` + +### Step 6 + +Only after the above succeeds, add the NRSS xarray metadata bridge. + +### Step 7 + +Promote the NRSS analytical comparison migration: + +- replace detector-space analytical tests with `NRSSIntegrator` / pure analytical checks where appropriate + +## Resume checklist for a fresh context + +If resuming from a new session, do this in order: + +1. read: + - `/homes/deand/dev/PyHyperScattering/docs/plans/nrss_q_coordinate_semantics.md` + - `/homes/deand/dev/NRSS/test-reports/nrss_q_coordinate_reduction_issue_draft.md` +2. confirm branch state: + - PyHyperScattering branch `221-unify-q-coordinate-semantics-between-nrss-simulation-outputs-and-pyhyperscattering-reduction-workflows` + - NRSS branch `backend/cupy-cyrsoxs-clone` +3. inspect current code anchors: + - `WPIntegrator.py` + - `cyrsoxsLoader.py` + - `morphology.py` + - `cupy_rsoxs.py` +4. implement `NRSSIntegrator` in PyHyperScattering first +5. add `test_NRSSIntegrator.py` +6. run PyHyperScattering tests +7. switch to `nrss-pyhyper-dev` +8. validate with local NRSS sphere outputs +9. add NRSS dev-only visualization tests +10. only then discuss NRSS metadata export changes + +## Commands likely to be useful + +These are not authoritative, just a restart aid. + +### PyHyperScattering tests + +```bash +pytest /homes/deand/dev/PyHyperScattering/tests/test_WPintegrator.py +pytest /homes/deand/dev/PyHyperScattering/tests/test_NRSSIntegrator.py +``` + +### NRSS analytical sphere reference + +```bash +pytest /homes/deand/dev/NRSS/tests/validation/test_analytical_sphere_form_factor.py -k sphere +``` + +### NRSS dev validation with plots + +```bash +NRSS_WRITE_VALIDATION_PLOTS=1 pytest /homes/deand/dev/NRSS/tests/validation/test_nrss_integrator_sphere_form_factor.py +``` + +## Non-goals for the first implementation + +- rewriting `WPIntegrator` +- replacing PyHyperScattering reduction with synthetic pyFAI geometry +- requiring NRSS metadata contract changes before testing +- migrating all NRSS validation files at once +- touching unrelated dirty NRSS worktree changes + +## Success criteria + +The work is successful when all of the following are true: + +1. PyHyperScattering exposes a working `NRSSIntegrator`. +2. For 2D NRSS output, `NRSSIntegrator` reproduces `WPIntegrator` semantics. +3. For 3D NRSS output, `NRSSIntegrator` emits corrected radial `q` semantics. +4. Local NRSS sphere validation shows `NRSSIntegrator` aligns with pure analytical form factor better than `WPIntegrator`. +5. Dev-only NRSS plots clearly show the semantic difference and support eventual migration of maintained tests. + diff --git a/docs/source/getting_started/integration.rst b/docs/source/getting_started/integration.rst index 2b15f30c..13d91620 100644 --- a/docs/source/getting_started/integration.rst +++ b/docs/source/getting_started/integration.rst @@ -23,7 +23,14 @@ For grazing-incidence experiments ``PGGeneralIntegrator`` uses ``pygix`` to output either reciprocal-space ``q_xy``/``q_z`` data or caked ``q`` vs ``chi``. Datasets already in ``qx``/``qy`` coordinates can be integrated with ``WPIntegrator`` which relies on -``skimage.transform.warp_polar`` (or its GPU version). +``skimage.transform.warp_polar`` (or its GPU version). ``WPIntegrator`` +is a detector-plane ``qx``/``qy`` to ``chi``/``q_perp`` remesher, so it +is the right choice for reciprocal-plane / 2D NRSS outputs. + +For NRSS detector-aware 3D outputs use ``NRSSIntegrator`` instead. +It keeps the same polar remesh approach but relabels the radial axis +with detector-corrected ``|q|`` so the reduced curve can be compared to +geometry-aware analytical form factors. All integrators return properly labeled xarray objects so the rest of the workflow can use normal xarray indexing and visualization. diff --git a/src/PyHyperScattering/Fitting.py b/src/PyHyperScattering/Fitting.py index 50dca688..5708c0db 100755 --- a/src/PyHyperScattering/Fitting.py +++ b/src/PyHyperScattering/Fitting.py @@ -94,15 +94,15 @@ def fit_lorentz(x,guess=None,pos_int_override=False,silent=False): except RuntimeError: if not silent: print("Fit failed to converge") - retval = xr.DataArray(data=np.nan,coords=x.coords).to_dataset(name='intensity') - retval['pos'] = xr.DataArray(data=np.nan,coords=x.coords) - retval['width'] = xr.DataArray(data=np.nan,coords=x.coords) + retval = xr.full_like(x, np.nan, dtype=float).to_dataset(name='intensity') + retval['pos'] = xr.full_like(x, np.nan, dtype=float) + retval['width'] = xr.full_like(x, np.nan, dtype=float) return retval if not silent: print(f"Fit completed, coeff = {coeff}") - retval = xr.DataArray(data=coeff[0],coords=x.coords).to_dataset(name='intensity') - retval['pos'] = xr.DataArray(data=coeff[1],coords=x.coords) - retval['width'] = xr.DataArray(data=coeff[2],coords=x.coords) + retval = xr.full_like(x, coeff[0], dtype=float).to_dataset(name='intensity') + retval['pos'] = xr.full_like(x, coeff[1], dtype=float) + retval['width'] = xr.full_like(x, coeff[2], dtype=float) return retval def fit_lorentz_bg(x,guess=None,pos_int_override=False,silent=False): ''' @@ -131,17 +131,17 @@ def fit_lorentz_bg(x,guess=None,pos_int_override=False,silent=False): except RuntimeError: if not silent: print("Fit failed to converge") - retval = xr.DataArray(data=np.nan,coords=x.coords).to_dataset(name='intensity') - retval['pos'] = xr.DataArray(data=np.nan,coords=x.coords) - retval['width'] = xr.DataArray(data=np.nan,coords=x.coords) - retval['bg'] = xr.DataArray(data=np.nan,coords=x.coords) + retval = xr.full_like(x, np.nan, dtype=float).to_dataset(name='intensity') + retval['pos'] = xr.full_like(x, np.nan, dtype=float) + retval['width'] = xr.full_like(x, np.nan, dtype=float) + retval['bg'] = xr.full_like(x, np.nan, dtype=float) return retval if not silent: print(f"Fit completed, coeff = {coeff}") - retval = xr.DataArray(data=coeff[0],coords=x.coords).to_dataset(name='intensity') - retval['pos'] = xr.DataArray(data=coeff[1],coords=x.coords) - retval['width'] = xr.DataArray(data=coeff[2],coords=x.coords) - retval['bg'] = xr.DataArray(data=coeff[3],coords=x.coords) + retval = xr.full_like(x, coeff[0], dtype=float).to_dataset(name='intensity') + retval['pos'] = xr.full_like(x, coeff[1], dtype=float) + retval['width'] = xr.full_like(x, coeff[2], dtype=float) + retval['bg'] = xr.full_like(x, coeff[3], dtype=float) return retval def fit_cos_anisotropy(data,qL,qU,qspacing,Enlist,ChiL,ChiU,binnumber,Chilim): @@ -329,4 +329,4 @@ def fit_cos(x_data, y_data): Ani_unc=0 Chisq=100 - return params, Ani, Ani_unc, Chisq \ No newline at end of file + return params, Ani, Ani_unc, Chisq diff --git a/src/PyHyperScattering/NRSSIntegrator.py b/src/PyHyperScattering/NRSSIntegrator.py new file mode 100644 index 00000000..26ffe4a1 --- /dev/null +++ b/src/PyHyperScattering/NRSSIntegrator.py @@ -0,0 +1,571 @@ +import ast +import warnings + +import numpy as np +import xarray as xr + +from PyHyperScattering.WPIntegrator import WPIntegrator + + +class NRSSIntegrator(WPIntegrator): + """ + Integrator for NRSS/CyRSoXS-style detector outputs stored on qx/qy coordinates. + + In reciprocal-plane / 2D mode the reduction is equivalent to WPIntegrator and the + radial coordinate is detector-plane q_perp. In detector-aware / 3D mode the same + polar remesh is used, but the radial coordinate is relabeled with detector-corrected + |q| based on the NRSS backend geometry. + """ + + _TWO_D_MODES = { + "2d", + "2d_reciprocal_plane", + "reciprocal_plane", + "reciprocal-plane", + "q_perp", + } + _THREE_D_MODES = { + "3d", + "3d_detector_aware", + "detector_aware", + "detector-aware", + "q_abs_detector_corrected", + } + + def __init__( + self, + return_cupy=False, + force_np_backend=False, + use_chunked_processing=False, + phys_size_nm=None, + shape_zyx=None, + z_dim=None, + energy_ev=None, + projection_mode=None, + validate_q_coords_against_phys_size="raise", + ): + super().__init__( + return_cupy=return_cupy, + force_np_backend=force_np_backend, + use_chunked_processing=use_chunked_processing, + ) + self._default_metadata = { + "phys_size_nm": phys_size_nm, + "shape_zyx": shape_zyx, + "z_dim": z_dim, + "energy_ev": energy_ev, + "projection_mode": projection_mode, + "validate_q_coords_against_phys_size": validate_q_coords_against_phys_size, + } + + def integrateSingleImage(self, img, **metadata_kwargs): + result = self.integrateImageStack_batched(img, **metadata_kwargs) + squeeze_dims = [dim for dim in result.dims if dim not in {"chi", "q"} and result.sizes[dim] == 1] + if squeeze_dims: + result = result.squeeze(squeeze_dims, drop=True) + return result + + def integrateImageStack(self, img_stack, method=None, chunksize=None, **metadata_kwargs): + if (self.use_chunked_processing and method is None) or method == "dask": + raise NotImplementedError( + "NRSSIntegrator does not support dask-backed reduction yet because " + "detector-corrected q coordinates can vary between slices." + ) + if method is None or method == "legacy" or method == "batched": + return self.integrateImageStack_batched(img_stack, **metadata_kwargs) + raise NotImplementedError(f"unsupported integration method {method}") + + def integrateImageStack_legacy(self, data, **metadata_kwargs): + return self.integrateImageStack_batched(data, **metadata_kwargs) + + def integrateImageStack_batched(self, data, **metadata_kwargs): + stacked, stacked_name, index_dims, spatial_dims = self._prepare_batched_input(data) + values = np.asarray(stacked.values) + if values.ndim != 3: + raise ValueError(f"NRSSIntegrator expected stacked data with 3 dims, got shape {values.shape!r}.") + + metadata = [self._resolve_metadata(stacked.isel({stacked_name: i}, drop=False), metadata_kwargs) for i in range(values.shape[0])] + reduced_values, q_axes, q_perp_axis, chi, common_attrs, common_mode, q_semantics_vary = ( + self._integrate_batched_array(stacked, values, metadata, spatial_dims) + ) + + result = xr.DataArray( + reduced_values, + dims=[stacked_name, "chi", "q"], + coords={ + stacked_name: stacked.coords[stacked_name], + "chi": chi, + }, + attrs=common_attrs, + ) + + result = self._assign_output_q_coords(result, stacked_name, q_axes, q_perp_axis, common_mode, q_semantics_vary) + + if len(index_dims) == 0: + return result.isel({stacked_name: 0}, drop=True) + if len(index_dims) > 1: + result = result.unstack(stacked_name) + result = result.transpose(*index_dims, "chi", "q") + return result + + + def _resolve_metadata(self, img, metadata_kwargs): + fallback = dict(self._default_metadata) + fallback.update(metadata_kwargs) + + shape_zyx = self._shape_from_attrs(img.attrs) + if shape_zyx is None: + shape_zyx = self._shape_from_value(fallback.get("shape_zyx")) + + z_dim = self._z_dim_from_attrs(img.attrs) + if z_dim is None and shape_zyx is not None: + z_dim = int(shape_zyx[0]) + if z_dim is None and fallback.get("z_dim") is not None: + z_dim = int(fallback["z_dim"]) + if z_dim is None and fallback.get("shape_zyx") is not None: + z_dim = int(self._shape_from_value(fallback["shape_zyx"])[0]) + + projection_mode = self._projection_mode_from_attrs(img.attrs) + if projection_mode is None and fallback.get("projection_mode") is not None: + projection_mode = str(fallback["projection_mode"]).strip().lower() + + phys_size_nm = self._phys_size_from_attrs(img.attrs) + if phys_size_nm is None and fallback.get("phys_size_nm") is not None: + phys_size_nm = float(fallback["phys_size_nm"]) + + validate_q_coords_against_phys_size = fallback.get("validate_q_coords_against_phys_size", "raise") + if phys_size_nm is not None and validate_q_coords_against_phys_size: + self._validate_q_axes_against_phys_size( + img=img, + phys_size_nm=phys_size_nm, + mode=validate_q_coords_against_phys_size, + ) + + energy_ev = self._energy_from_attrs_or_coords(img) + if energy_ev is None and fallback.get("energy_ev") is not None: + energy_ev = float(fallback["energy_ev"]) + + nrss_semantic_mode = self._semantic_mode(projection_mode=projection_mode, z_dim=z_dim) + radial_semantics = "q_perp" + if nrss_semantic_mode == "3d_detector_aware": + radial_semantics = "q_abs_detector_corrected" + if energy_ev is None: + raise ValueError( + "NRSSIntegrator needs a scalar energy for detector-aware 3D reductions. " + "Populate img.attrs['energy_ev'], provide a scalar energy coordinate, or pass energy_ev=..." + ) + + return { + "phys_size_nm": phys_size_nm, + "shape_zyx": shape_zyx, + "z_dim": z_dim, + "energy_ev": energy_ev, + "nrss_semantic_mode": nrss_semantic_mode, + "radial_semantics": radial_semantics, + } + + @staticmethod + def _spatial_dims(img): + missing = [name for name in ("qx", "qy") if name not in img.coords] + if missing: + raise ValueError(f"NRSSIntegrator requires qx/qy coordinates, missing {missing}.") + dims = tuple(dim for dim in img.dims if dim in {"qx", "qy"}) + if set(dims) != {"qx", "qy"}: + warnings.warn( + "NRSSIntegrator found qx/qy coordinates but not both as dimensions. " + "Reduction will proceed using the available qx/qy coordinates.", + stacklevel=2, + ) + return ("qx", "qy") + + @staticmethod + def _axis_center(axis, name): + return float( + xr.DataArray(np.linspace(0, len(axis) - 1, len(axis))) + .assign_coords({"dim_0": axis.values}) + .rename({"dim_0": name}) + .interp({name: 0}) + .data + ) + + @staticmethod + def _q_perp_axis(img, n_points): + q = np.sqrt(img.qy ** 2 + img.qx ** 2) + return np.linspace(0.0, float(np.nanmax(q)), int(n_points), dtype=np.float64) + + @staticmethod + def _detector_corrected_q(q_perp_axis, energy_ev): + q_perp_axis = np.asarray(q_perp_axis, dtype=np.float64) + wavelength_nm = 1239.84197 / float(energy_ev) + k = 2.0 * np.pi / wavelength_nm + val = k * k - q_perp_axis * q_perp_axis + qz = np.full_like(q_perp_axis, np.nan, dtype=np.float64) + q = np.full_like(q_perp_axis, np.nan, dtype=np.float64) + valid = val >= 0.0 + qz[valid] = -k + np.sqrt(val[valid]) + q[valid] = np.sqrt(q_perp_axis[valid] * q_perp_axis[valid] + qz[valid] * qz[valid]) + return q + + @staticmethod + def _expected_detector_axis(n_points, phys_size_nm): + return 2.0 * np.pi * np.fft.fftshift(np.fft.fftfreq(int(n_points), d=float(phys_size_nm))) + + def _validate_q_axes_against_phys_size(self, img, phys_size_nm, mode): + qx = np.asarray(img.qx.values, dtype=np.float64) + qy = np.asarray(img.qy.values, dtype=np.float64) + expected_qx = self._expected_detector_axis(qx.size, phys_size_nm) + expected_qy = self._expected_detector_axis(qy.size, phys_size_nm) + qx_ok = np.allclose(qx, expected_qx, atol=1e-12, rtol=0.0, equal_nan=True) + qy_ok = np.allclose(qy, expected_qy, atol=1e-12, rtol=0.0, equal_nan=True) + if qx_ok and qy_ok: + return + + message = ( + "NRSSIntegrator detected inconsistent q coordinates for the provided phys_size_nm. " + "qx/qy are authoritative for reduction, but they do not match the FFT-style detector " + "axes implied by phys_size_nm." + ) + if mode == "raise": + raise ValueError(message) + if mode == "warn": + warnings.warn(message, stacklevel=2) + return + raise ValueError( + "validate_q_coords_against_phys_size must be one of False, 'warn', or 'raise'. " + f"Got {mode!r}." + ) + + def _semantic_mode(self, projection_mode, z_dim): + if projection_mode is not None: + mode = str(projection_mode).strip().lower() + if mode in self._TWO_D_MODES: + return "2d_reciprocal_plane" + if mode in self._THREE_D_MODES: + return "3d_detector_aware" + raise ValueError(f"Unrecognized NRSS projection_mode {projection_mode!r}.") + if z_dim is None: + raise ValueError( + "NRSSIntegrator could not determine whether the input is 2D reciprocal-plane or " + "3D detector-aware. Provide z_dim, shape_zyx, or projection_mode metadata." + ) + return "2d_reciprocal_plane" if int(z_dim) == 1 else "3d_detector_aware" + + @staticmethod + def _shape_from_attrs(attrs): + for key in ("shape_zyx", "num_zyx", "NumZYX"): + if key in attrs: + return NRSSIntegrator._shape_from_value(attrs[key]) + return None + + @staticmethod + def _z_dim_from_attrs(attrs): + if "z_dim" in attrs: + return int(attrs["z_dim"]) + shape_zyx = NRSSIntegrator._shape_from_attrs(attrs) + if shape_zyx is None: + return None + return int(shape_zyx[0]) + + @staticmethod + def _phys_size_from_attrs(attrs): + for key in ("phys_size_nm", "PhysSize"): + if key in attrs and attrs[key] is not None: + return float(attrs[key]) + return None + + @staticmethod + def _projection_mode_from_attrs(attrs): + for key in ("nrss_output_semantics", "nrss_semantic_mode", "projection_mode"): + value = attrs.get(key) + if value is not None: + return str(value).strip().lower() + return None + + @staticmethod + def _energy_from_attrs_or_coords(img): + for key in ("energy_ev", "Energy"): + value = img.attrs.get(key) + if value is not None: + return float(value) + + if "energy" not in img.coords: + return None + + energy_values = np.asarray(img.coords["energy"].values, dtype=np.float64).reshape(-1) + if energy_values.size == 1: + return float(energy_values[0]) + return None + + @staticmethod + def _shape_from_value(value): + if value is None: + return None + if isinstance(value, str): + value = ast.literal_eval(value) + shape = tuple(int(v) for v in value) + if len(shape) != 3: + raise ValueError(f"shape_zyx/num_zyx must contain exactly 3 entries, got {shape!r}.") + return shape + + @staticmethod + def _prepare_batched_input(data): + spatial_dims = NRSSIntegrator._spatial_dims(data) + index_dims = [dim for dim in data.dims if dim not in spatial_dims] + spatial_dims_in_order = tuple(dim for dim in data.dims if dim in {"qx", "qy"}) + + if len(index_dims) == 0: + stacked_name = "pyhyper_internal_batch" + stacked = data.expand_dims({stacked_name: [0]}) + elif len(index_dims) == 1: + stacked_name = index_dims[0] + stacked = data + else: + stacked_name = "pyhyper_internal_multiindex" + stacked = data.stack({stacked_name: index_dims}) + + stacked = stacked.transpose(stacked_name, *spatial_dims_in_order) + return stacked, stacked_name, index_dims, spatial_dims_in_order + + def _integrate_batched_array(self, stacked, values, metadata, spatial_dims): + center_x = self._axis_center(stacked.qx, "qx") + center_y = self._axis_center(stacked.qy, "qy") + center_lookup = { + "qx": center_x, + "qy": center_y, + } + center = tuple(center_lookup[dim] for dim in spatial_dims) + radius = np.sqrt((values.shape[1] - center[0]) ** 2 + (values.shape[2] - center[1]) ** 2) + reduced_values = self._warp_polar_batched(values, center=center, radius=radius) + + q_perp_axis = self._q_perp_axis(stacked, int(reduced_values.shape[-1])) + q_axes, common_mode, q_semantics_vary = self._batched_q_axes(q_perp_axis, metadata) + chi = np.linspace(-179.5, 179.5, reduced_values.shape[1]) + common_attrs = self._build_batched_attrs(stacked.attrs, metadata[0]) + return reduced_values, q_axes, q_perp_axis, chi, common_attrs, common_mode, q_semantics_vary + + def _assign_output_q_coords(self, result, stacked_name, q_axes, q_perp_axis, common_mode, q_semantics_vary): + q_axes_are_same = self._allclose_1d(q_axes) + if q_axes_are_same and not q_semantics_vary: + result = result.assign_coords(q=q_axes[0]) + if q_perp_axis is not None and not np.allclose(q_axes[0], q_perp_axis, atol=0.0, rtol=0.0): + result = result.assign_coords(q_perp=("q", q_perp_axis)) + return result + + if common_mode == "3d_detector_aware" and not q_semantics_vary: + q_common = self._shared_q_grid(q_axes) + if q_common is not None: + interpolated = self._interp_stack_to_common_q(np.asarray(result.values), q_axes, q_common) + result = xr.DataArray( + interpolated, + dims=result.dims, + coords={ + stacked_name: result.coords[stacked_name], + "chi": result.coords["chi"].values, + "q": q_common, + "q_abs": ((stacked_name, "q"), np.stack(q_axes, axis=0)), + }, + attrs=dict(result.attrs), + ) + result.attrs["radial_coordinate_mode"] = "shared_q_grid_interpolated" + result.attrs["q_axis_note"] = ( + "The q dimension is a shared detector-corrected q grid spanning the overlap " + "of all slices. Exact per-slice q values before interpolation remain in q_abs." + ) + result.attrs.pop("energy_ev", None) + return result + + return self._attach_per_slice_q_index(result, stacked_name, q_axes, q_perp_axis) + + @staticmethod + def _build_batched_attrs(base_attrs, metadata): + attrs = dict(base_attrs) + attrs.update( + { + "radial_semantics": metadata["radial_semantics"], + "source_integrator": "NRSSIntegrator", + "nrss_semantic_mode": metadata["nrss_semantic_mode"], + "phys_size_nm": metadata["phys_size_nm"], + "z_dim": metadata["z_dim"], + } + ) + if metadata["shape_zyx"] is not None: + attrs["shape_zyx"] = tuple(metadata["shape_zyx"]) + if metadata["energy_ev"] is not None: + attrs["energy_ev"] = float(metadata["energy_ev"]) + return attrs + + @staticmethod + def _batched_q_axes(q_perp_axis, metadata): + q_axes = [] + common_mode = None + q_semantics_vary = False + three_d_indices = [] + three_d_energies = [] + + for i, md in enumerate(metadata): + mode = md["nrss_semantic_mode"] + if common_mode is None: + common_mode = mode + elif mode != common_mode: + q_semantics_vary = True + q_axes.append(np.asarray(q_perp_axis, dtype=np.float64)) + if mode == "3d_detector_aware": + three_d_indices.append(i) + three_d_energies.append(md["energy_ev"]) + + if three_d_indices: + corrected = NRSSIntegrator._detector_corrected_q_batch(q_perp_axis, np.asarray(three_d_energies, dtype=np.float64)) + for idx, q_axis in zip(three_d_indices, corrected): + q_axes[idx] = q_axis + + return q_axes, common_mode, q_semantics_vary + + def _warp_polar_batched(self, values, center, radius): + if self.MACHINE_HAS_CUDA: + try: + import cupy as cp + except ImportError: # pragma: no cover + return self._warp_polar_batched_numpy(values, center=center, radius=radius) + + values_xp = cp.asarray(values) + reduced = self._warp_polar_batched_xp(values_xp, center=center, radius=radius, xp=cp) + if self.return_cupy: + return reduced + return cp.asnumpy(reduced) + + return self._warp_polar_batched_numpy(values, center=center, radius=radius) + + def _warp_polar_batched_numpy(self, values, center, radius): + reduced = self._warp_polar_batched_xp(np.asarray(values), center=center, radius=radius, xp=np) + return np.asarray(reduced) + + @staticmethod + def _warp_polar_batched_xp(values, center, radius, xp): + values = xp.asarray(values) + n_images, n_rows, n_cols = values.shape + n_theta = 360 + n_radius = int(np.ceil(radius)) + if n_radius <= 0: + raise ValueError(f"NRSSIntegrator computed a non-positive polar radius {radius!r}.") + + center_row, center_col = center + theta = xp.deg2rad(xp.arange(n_theta, dtype=xp.float64)) + radial = xp.arange(n_radius, dtype=xp.float64) * (float(radius) / n_radius) + radial_grid, theta_grid = xp.meshgrid(radial, theta) + + row_coords = radial_grid * xp.sin(theta_grid) + center_row + col_coords = radial_grid * xp.cos(theta_grid) + center_col + + row0 = xp.floor(row_coords).astype(xp.int64) + col0 = xp.floor(col_coords).astype(xp.int64) + row1 = row0 + 1 + col1 = col0 + 1 + + row_weight = row_coords - row0 + col_weight = col_coords - col0 + + def sample(row_idx, col_idx): + valid = (row_idx >= 0) & (row_idx < n_rows) & (col_idx >= 0) & (col_idx < n_cols) + row_clip = xp.clip(row_idx, 0, n_rows - 1) + col_clip = xp.clip(col_idx, 0, n_cols - 1) + sampled = values[:, row_clip, col_clip] + return sampled * valid[None, :, :] + + top_left = sample(row0, col0) + top_right = sample(row0, col1) + bottom_left = sample(row1, col0) + bottom_right = sample(row1, col1) + + return ( + top_left * (1.0 - row_weight)[None, :, :] * (1.0 - col_weight)[None, :, :] + + top_right * (1.0 - row_weight)[None, :, :] * col_weight[None, :, :] + + bottom_left * row_weight[None, :, :] * (1.0 - col_weight)[None, :, :] + + bottom_right * row_weight[None, :, :] * col_weight[None, :, :] + ) + + @staticmethod + def _detector_corrected_q_batch(q_perp_axis, energy_ev): + q_perp_axis = np.asarray(q_perp_axis, dtype=np.float64)[None, :] + energy_ev = np.asarray(energy_ev, dtype=np.float64).reshape(-1, 1) + wavelength_nm = 1239.84197 / energy_ev + k = 2.0 * np.pi / wavelength_nm + val = k * k - q_perp_axis * q_perp_axis + valid = val >= 0.0 + qz = -k + np.sqrt(val, where=valid, out=np.full_like(val, np.nan, dtype=np.float64)) + q = np.full_like(val, np.nan, dtype=np.float64) + q_perp_broadcast = np.broadcast_to(q_perp_axis, val.shape) + q[valid] = np.sqrt(q_perp_broadcast[valid] * q_perp_broadcast[valid] + qz[valid] * qz[valid]) + return q + + @staticmethod + def _allclose_1d(arrays): + if len(arrays) <= 1: + return True + ref = np.asarray(arrays[0], dtype=np.float64) + for arr in arrays[1:]: + candidate = np.asarray(arr, dtype=np.float64) + if candidate.shape != ref.shape: + return False + if not np.allclose(candidate, ref, equal_nan=True, rtol=0.0, atol=1e-12): + return False + return True + + @staticmethod + def _shared_q_grid(q_axes): + lower_bounds = [] + upper_bounds = [] + n_points = None + + for axis in q_axes: + axis = np.asarray(axis, dtype=np.float64) + finite = axis[np.isfinite(axis)] + if finite.size < 2: + return None + lower_bounds.append(float(np.min(finite))) + upper_bounds.append(float(np.max(finite))) + n_points = axis.size if n_points is None else min(n_points, axis.size) + + q_min = max(lower_bounds) + q_max = min(upper_bounds) + if not np.isfinite(q_min) or not np.isfinite(q_max) or q_max <= q_min or n_points is None or n_points < 2: + return None + return np.linspace(q_min, q_max, int(n_points), dtype=np.float64) + + @staticmethod + def _interp_stack_to_common_q(values, q_axes, q_common): + values = np.asarray(values) + output = np.full((values.shape[0], values.shape[1], q_common.size), np.nan, dtype=values.dtype) + + for i, q_axis in enumerate(q_axes): + q_axis = np.asarray(q_axis, dtype=np.float64) + valid = np.isfinite(q_axis) + q_valid = q_axis[valid] + if q_valid.size < 2: + continue + + q_valid, unique_idx = np.unique(q_valid, return_index=True) + slice_values = values[i][:, valid][:, unique_idx] + for chi_idx in range(values.shape[1]): + output[i, chi_idx, :] = np.interp( + q_common, + q_valid, + slice_values[chi_idx], + left=np.nan, + right=np.nan, + ) + + return output + + @staticmethod + def _attach_per_slice_q_index(result, stacked_name, q_axes, q_perp_axis): + result = result.assign_coords(q=np.arange(result.sizes["q"], dtype=np.int64)) + result = result.assign_coords(q_abs=((stacked_name, "q"), np.stack(q_axes, axis=0))) + if q_perp_axis is not None: + result = result.assign_coords(q_perp=("q", q_perp_axis)) + result.attrs["radial_coordinate_mode"] = "per_slice_q_abs" + result.attrs["q_axis_note"] = ( + "The q dimension indexes radial bins. Exact detector-corrected q values are " + "stored in the q_abs coordinate." + ) + result.attrs.pop("energy_ev", None) + return result diff --git a/src/PyHyperScattering/PFGeneralIntegrator.py b/src/PyHyperScattering/PFGeneralIntegrator.py index 26a95c08..b5844c2a 100755 --- a/src/PyHyperScattering/PFGeneralIntegrator.py +++ b/src/PyHyperScattering/PFGeneralIntegrator.py @@ -1,6 +1,7 @@ from pyFAI import azimuthalIntegrator from pyFAI.units import eq_q, formula_q, register_radial_unit from pyFAI.io.ponifile import PoniFile +import io import h5py import warnings import xarray as xr @@ -555,7 +556,10 @@ def loadPyHyperMask(self, **kwargs): # print(strlist) dflist = [] for item in strlist: - dflist.append(pd.read_json(item)) + if isinstance(item, str): + dflist.append(pd.read_json(io.StringIO(item))) + else: + dflist.append(pd.DataFrame(item)) # print(dflist) pyhyperlist = [] for shape in dflist: diff --git a/src/PyHyperScattering/SST1RSoXSDB.py b/src/PyHyperScattering/SST1RSoXSDB.py index 5148ccb4..56c30db5 100644 --- a/src/PyHyperScattering/SST1RSoXSDB.py +++ b/src/PyHyperScattering/SST1RSoXSDB.py @@ -32,6 +32,15 @@ import copy +def _scalarize_dark_index(value): + arr = np.asarray(value).reshape(-1) + if arr.size == 0: + raise ValueError("dark_id is empty") + if not np.all(arr == arr[0]): + raise ValueError(f"dark_id is not uniform within group: {arr}") + return int(arr[0]) + + class SST1RSoXSDB: """ Loader for bluesky run xarrays form NSLS-II SST1 RSoXS instrument @@ -906,7 +915,8 @@ def loadRun( data = data.assign_coords(dark_id=("time", darkframe)) def subtract_dark(img, pedestal=100, darks=None): - return img + pedestal - darks[int(img.dark_id.values)] + dark_idx = _scalarize_dark_index(img.dark_id.values) + return img + pedestal - darks[dark_idx] data = data.groupby("time",squeeze=False).map(subtract_dark, darks=dark, pedestal=self.dark_pedestal) diff --git a/src/PyHyperScattering/WPIntegrator.py b/src/PyHyperScattering/WPIntegrator.py index 8b4ed981..01483509 100644 --- a/src/PyHyperScattering/WPIntegrator.py +++ b/src/PyHyperScattering/WPIntegrator.py @@ -26,7 +26,12 @@ class WPIntegrator(): ''' - Integrator for qx/qy format xarrays using skimage.transform.warp_polar or a custom cuda-accelerated version, warp_polar_gpu + Integrator for qx/qy format xarrays using skimage.transform.warp_polar or a custom + cuda-accelerated version, warp_polar_gpu. + + The radial output of this integrator is detector-plane q_perp. This is the correct + semantics for reciprocal-plane / 2D NRSS outputs, but not for detector-aware 3D + NRSS outputs that should be compared against detector-corrected |q|. ''' def __init__(self,return_cupy=False,force_np_backend=False,use_chunked_processing=False): @@ -101,15 +106,17 @@ def integrateSingleImage(self,img): .rename({'dim_0':'qy'}) .interp(qy=0) .data) + stacked_axis = None + system_to_integ = None + candidate_axes = list(img.coords) try: - stacked_axis = list(img.coords) - stacked_axis.remove('qx') - stacked_axis.remove('qy') - assert len(stacked_axis)==1, f"More than one axis left ({stacked_axis}) after removing qx and qy, not sure how to handle" - stacked_axis = stacked_axis[0] + candidate_axes.remove('qx') + candidate_axes.remove('qy') + except ValueError: + candidate_axes = [] + if len(candidate_axes) == 1: + stacked_axis = candidate_axes[0] system_to_integ = img.__getattr__(stacked_axis) - except AttributeError: - pass if self.MACHINE_HAS_CUDA: TwoD = self.warp_polar_gpu(img_to_integ,center=(center_x,center_y), radius = np.sqrt((img_to_integ.shape[0] - center_x)**2 + (img_to_integ.shape[1] - center_y)**2)) @@ -126,10 +133,19 @@ def integrateSingleImage(self,img): chi = np.linspace(-179.5,179.5,360) # chi = np.linspace(0.5,359.5,360) + attrs = dict(img.attrs) + attrs.update( + { + "radial_semantics": "q_perp", + "source_integrator": "WPIntegrator", + } + ) try: - return xr.DataArray([TwoD],dims=[stacked_axis,'chi','q'],coords={'q':q,'chi':chi,stacked_axis:system_to_integ},attrs=img.attrs) + if stacked_axis is None or system_to_integ is None: + raise ValueError("No stacked axis present.") + return xr.DataArray([TwoD],dims=[stacked_axis,'chi','q'],coords={'q':q,'chi':chi,stacked_axis:system_to_integ},attrs=attrs) except ValueError: - return xr.DataArray(TwoD,dims=['chi','q'],coords={'q':q,'chi':chi},attrs=img.attrs) + return xr.DataArray(TwoD,dims=['chi','q'],coords={'q':q,'chi':chi},attrs=attrs) def integrateImageStack(self,img_stack,method=None,chunksize=None): @@ -160,6 +176,9 @@ def integrateImageStack_legacy(self,data): indexes.remove('qy') except ValueError: pass + + if len(indexes) == 0: + return self.integrateSingleImage(data) if len(indexes) == 1: if data.__getattr__(indexes[0]).to_pandas().drop_duplicates().shape[0] != data.__getattr__(indexes[0]).shape[0]: @@ -189,6 +208,9 @@ def integrateImageStack_dask(self,data,chunksize=5): indexes.remove('qy') except ValueError: pass + + if len(indexes) == 0: + return self.integrateSingleImage(data) if len(indexes) == 1: if data.__getattr__(indexes[0]).to_pandas().drop_duplicates().shape[0] != data.__getattr__(indexes[0]).shape[0]: diff --git a/src/PyHyperScattering/integrate.py b/src/PyHyperScattering/integrate.py index 94000a81..193ae9cf 100644 --- a/src/PyHyperScattering/integrate.py +++ b/src/PyHyperScattering/integrate.py @@ -1,4 +1,5 @@ from PyHyperScattering.PFEnergySeriesIntegrator import PFEnergySeriesIntegrator from PyHyperScattering.PFGeneralIntegrator import PFGeneralIntegrator from PyHyperScattering.WPIntegrator import WPIntegrator -from PyHyperScattering.PGGeneralIntegrator import PGGeneralIntegrator \ No newline at end of file +from PyHyperScattering.PGGeneralIntegrator import PGGeneralIntegrator +from PyHyperScattering.NRSSIntegrator import NRSSIntegrator diff --git a/tests/test_Fitting_lorentz.py b/tests/test_Fitting_lorentz.py index 53f59ad2..d5ab3904 100644 --- a/tests/test_Fitting_lorentz.py +++ b/tests/test_Fitting_lorentz.py @@ -26,6 +26,8 @@ def test_fit_lorentz(): res = Fitting.fit_lorentz( lorentz_da, guess=[AMPLITUDE * 0.9, CENTER * 1.01, WIDTH * 1.1], silent=True ) + assert list(res.intensity.dims) == ["q"] + assert np.allclose(res.intensity.values, AMPLITUDE, rtol=1e-5) assert np.isclose(res.intensity.mean(), AMPLITUDE, rtol=1e-5) assert np.isclose(res.pos.mean(), CENTER, rtol=1e-5) assert np.isclose(res.width.mean(), WIDTH, rtol=1e-5) @@ -53,3 +55,19 @@ def test_fitting_apply_with_stacked_dimension(): assert list(result.dims) == ["replicate"] assert list(result.coords["replicate"].values) == ["r1", "r2"] +def test_fit_lorentz_output_is_plottable(): + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + res = Fitting.fit_lorentz_bg( + lorentz_bg_da, + guess=[AMPLITUDE * 0.9, CENTER * 1.01, WIDTH * 1.1, BACKGROUND * 0.8], + silent=True, + ) + + fig, ax = plt.subplots() + artist = res.intensity.plot(ax=ax) + assert artist is not None + assert list(res.intensity.dims) == ["q"] + plt.close(fig) diff --git a/tests/test_NRSSIntegrator.py b/tests/test_NRSSIntegrator.py new file mode 100644 index 00000000..414c406d --- /dev/null +++ b/tests/test_NRSSIntegrator.py @@ -0,0 +1,213 @@ +import numpy as np +import pytest +import xarray as xr + +from PyHyperScattering.integrate import NRSSIntegrator, WPIntegrator + + +def _synthetic_qxy_image( + *, + nx=101, + ny=91, + qmax=1.0, + phys_size_nm=None, + attrs=None, + dims=("qy", "qx"), +): + if phys_size_nm is None: + qx = np.linspace(-qmax, qmax, nx, dtype=np.float64) + qy = np.linspace(-qmax, qmax, ny, dtype=np.float64) + else: + qx = 2.0 * np.pi * np.fft.fftshift(np.fft.fftfreq(nx, d=float(phys_size_nm))) + qy = 2.0 * np.pi * np.fft.fftshift(np.fft.fftfreq(ny, d=float(phys_size_nm))) + qx_grid, qy_grid = np.meshgrid(qx, qy) + image = ( + np.exp(-((qx_grid - 0.18) ** 2 + (qy_grid + 0.07) ** 2) / 0.010) + + 0.65 * np.exp(-((qx_grid + 0.24) ** 2 + (qy_grid - 0.22) ** 2) / 0.018) + + 0.08 + ) + da = xr.DataArray( + image, + dims=("qy", "qx"), + coords={"qy": qy, "qx": qx}, + attrs=dict(attrs or {}), + ) + if tuple(dims) == ("qy", "qx"): + return da + if tuple(dims) == ("qx", "qy"): + return da.transpose("qx", "qy") + raise ValueError(f"Unsupported dims ordering {dims!r}.") + + +def test_nrss_integrator_2d_matches_wp_semantics(): + img = _synthetic_qxy_image( + phys_size_nm=5.0, + attrs={"z_dim": 1, "phys_size_nm": 5.0}, + dims=("qx", "qy"), + ) + nrss = NRSSIntegrator(force_np_backend=True) + wp = WPIntegrator(force_np_backend=True) + + reduced_nrss = nrss.integrateImageStack(img) + reduced_wp = wp.integrateImageStack(img) + + np.testing.assert_allclose(reduced_nrss.coords["q"].values, reduced_wp.coords["q"].values) + np.testing.assert_allclose(reduced_nrss.values, reduced_wp.values) + assert reduced_nrss.attrs["radial_semantics"] == "q_perp" + assert reduced_nrss.attrs["nrss_semantic_mode"] == "2d_reciprocal_plane" + assert reduced_nrss.attrs["source_integrator"] == "NRSSIntegrator" + + +def test_nrss_integrator_3d_corrects_radial_axis(): + energy_ev = 285.0 + img = _synthetic_qxy_image( + phys_size_nm=5.0, + attrs={"z_dim": 128, "phys_size_nm": 5.0, "energy_ev": energy_ev}, + ) + reduced = NRSSIntegrator(force_np_backend=True).integrateImageStack(img) + + q = np.asarray(reduced.coords["q"].values, dtype=np.float64) + q_perp = np.asarray(reduced.coords["q_perp"].values, dtype=np.float64) + wavelength_nm = 1239.84197 / energy_ev + k = 2.0 * np.pi / wavelength_nm + expected = np.full_like(q_perp, np.nan) + valid = (k * k - q_perp * q_perp) >= 0.0 + qz = -k + np.sqrt(k * k - q_perp[valid] * q_perp[valid]) + expected[valid] = np.sqrt(q_perp[valid] * q_perp[valid] + qz * qz) + + np.testing.assert_allclose(q[valid], expected[valid], rtol=0.0, atol=1e-12) + assert np.nanmax(np.abs(q[valid] - q_perp[valid])) > 1e-3 + assert reduced.attrs["radial_semantics"] == "q_abs_detector_corrected" + assert reduced.attrs["nrss_semantic_mode"] == "3d_detector_aware" + + +def test_nrss_integrator_falls_back_to_explicit_kwargs(): + img = _synthetic_qxy_image(phys_size_nm=5.0, attrs={}) + reduced = NRSSIntegrator(force_np_backend=True).integrateImageStack( + img, + phys_size_nm=5.0, + z_dim=64, + energy_ev=285.0, + ) + + assert isinstance(reduced, xr.DataArray) + assert reduced.attrs["phys_size_nm"] == 5.0 + assert reduced.attrs["z_dim"] == 64 + assert reduced.attrs["energy_ev"] == 285.0 + assert reduced.attrs["nrss_semantic_mode"] == "3d_detector_aware" + + +def test_nrss_integrator_preserves_stack_axis(): + img0 = _synthetic_qxy_image( + phys_size_nm=5.0, + attrs={"z_dim": 1, "phys_size_nm": 5.0}, + dims=("qx", "qy"), + ) + img1 = (1.15 * img0).assign_coords(qx=img0.qx, qy=img0.qy) + stacked = xr.concat([img0, img1], dim=xr.IndexVariable("energy", [285.0, 286.0])) + + reduced = NRSSIntegrator(force_np_backend=True).integrateImageStack(stacked) + + assert reduced.dims == ("energy", "chi", "q") + assert list(reduced.coords["energy"].values) == [285.0, 286.0] + assert reduced.sizes["chi"] == 360 + assert reduced.attrs["radial_semantics"] == "q_perp" + + +def test_nrss_integrator_batched_matches_legacy_for_2d_stack(): + img0 = _synthetic_qxy_image( + phys_size_nm=5.0, + attrs={"z_dim": 1, "phys_size_nm": 5.0}, + dims=("qx", "qy"), + ) + img1 = (1.15 * img0).assign_coords(qx=img0.qx, qy=img0.qy) + stacked = xr.concat([img0, img1], dim=xr.IndexVariable("energy", [285.0, 286.0])) + + integrator = NRSSIntegrator(force_np_backend=True) + legacy = integrator.integrateImageStack(stacked, method="legacy") + batched = integrator.integrateImageStack(stacked, method="batched") + + np.testing.assert_allclose(batched.values, legacy.values, rtol=1e-12, atol=1e-12) + np.testing.assert_allclose(batched.q.values, legacy.q.values, rtol=0.0, atol=1e-12) + + +def test_nrss_integrator_3d_stack_uses_shared_physical_q_axis(): + img0 = _synthetic_qxy_image( + phys_size_nm=5.0, + attrs={"z_dim": 64, "phys_size_nm": 5.0}, + dims=("qx", "qy"), + ) + img1 = (1.15 * img0).assign_coords(qx=img0.qx, qy=img0.qy) + stacked = xr.concat([img0, img1], dim=xr.IndexVariable("energy", [285.0, 286.0])) + + reduced = NRSSIntegrator(force_np_backend=True).integrateImageStack(stacked) + + q = np.asarray(reduced.coords["q"].values, dtype=np.float64) + assert reduced.dims == ("energy", "chi", "q") + assert reduced.attrs["radial_semantics"] == "q_abs_detector_corrected" + assert reduced.attrs["radial_coordinate_mode"] == "shared_q_grid_interpolated" + assert "q_abs" in reduced.coords + assert np.issubdtype(q.dtype, np.floating) + assert np.all(np.diff(q) > 0.0) + assert not np.array_equal(q, np.arange(q.size, dtype=np.float64)) + + +def test_nrss_integrator_batched_matches_legacy_for_3d_stack(): + img0 = _synthetic_qxy_image( + phys_size_nm=5.0, + attrs={"z_dim": 64, "phys_size_nm": 5.0}, + dims=("qx", "qy"), + ) + img1 = (1.15 * img0).assign_coords(qx=img0.qx, qy=img0.qy) + stacked = xr.concat([img0, img1], dim=xr.IndexVariable("energy", [285.0, 286.0])) + + integrator = NRSSIntegrator(force_np_backend=True) + legacy = integrator.integrateImageStack(stacked, method="legacy") + batched = integrator.integrateImageStack(stacked, method="batched") + + np.testing.assert_allclose(batched.values, legacy.values, rtol=1e-7, atol=1e-7) + np.testing.assert_allclose(batched.q.values, legacy.q.values, rtol=0.0, atol=1e-12) + np.testing.assert_allclose(batched.q_abs.values, legacy.q_abs.values, rtol=0.0, atol=1e-12) + + +def test_wp_integrator_semantics_note(): + img = _synthetic_qxy_image() + reduced = WPIntegrator(force_np_backend=True).integrateImageStack(img) + + assert reduced.attrs["radial_semantics"] == "q_perp" + assert reduced.attrs["source_integrator"] == "WPIntegrator" + + +def test_wp_integrator_preserves_minimal_nrss_metadata_attrs(): + img = _synthetic_qxy_image(phys_size_nm=5.0, attrs={"phys_size_nm": 5.0, "z_dim": 64}) + reduced = WPIntegrator(force_np_backend=True).integrateImageStack(img) + + assert reduced.attrs["phys_size_nm"] == 5.0 + assert reduced.attrs["z_dim"] == 64 + assert reduced.attrs["radial_semantics"] == "q_perp" + + +def test_nrss_integrator_raises_on_phys_size_q_axis_mismatch_by_default(): + img = _synthetic_qxy_image(attrs={"phys_size_nm": 5.0, "z_dim": 64, "energy_ev": 285.0}) + with np.testing.assert_raises_regex( + ValueError, + "inconsistent q coordinates", + ): + NRSSIntegrator(force_np_backend=True).integrateImageStack(img) + + +def test_nrss_integrator_can_warn_or_skip_phys_size_q_axis_validation(): + img = _synthetic_qxy_image(attrs={"phys_size_nm": 5.0, "z_dim": 64, "energy_ev": 285.0}) + + with pytest.warns(UserWarning, match="inconsistent q coordinates"): + reduced_warn = NRSSIntegrator( + force_np_backend=True, + validate_q_coords_against_phys_size="warn", + ).integrateImageStack(img) + assert isinstance(reduced_warn, xr.DataArray) + + reduced_skip = NRSSIntegrator( + force_np_backend=True, + validate_q_coords_against_phys_size=False, + ).integrateImageStack(img) + assert isinstance(reduced_skip, xr.DataArray) diff --git a/tests/test_PFGeneralIntegrator_mask_parsing.py b/tests/test_PFGeneralIntegrator_mask_parsing.py new file mode 100644 index 00000000..821d3e80 --- /dev/null +++ b/tests/test_PFGeneralIntegrator_mask_parsing.py @@ -0,0 +1,24 @@ +import json +import pathlib + +import pandas as pd + +from PyHyperScattering.PFGeneralIntegrator import PFGeneralIntegrator + + +def test_load_pyhyper_mask_parses_json_strings_with_pandas3(tmp_path): + polygon_df = pd.DataFrame( + { + "x": [2, 7, 7, 2], + "y": [2, 2, 7, 7], + } + ) + mask_file = tmp_path / "mask.json" + mask_file.write_text(json.dumps([polygon_df.to_json()])) + + integrator = PFGeneralIntegrator(maskmethod="none", geomethod="none") + integrator.loadPyHyperMask(maskpath=pathlib.Path(mask_file), maskshape=(10, 10)) + + assert integrator.mask.shape == (10, 10) + assert integrator.mask.dtype == bool + assert integrator.mask.any() diff --git a/tests/test_SST1RSoXSDB_dark_index.py b/tests/test_SST1RSoXSDB_dark_index.py new file mode 100644 index 00000000..3092ad36 --- /dev/null +++ b/tests/test_SST1RSoXSDB_dark_index.py @@ -0,0 +1,21 @@ +import numpy as np +import pytest + +from PyHyperScattering.SST1RSoXSDB import _scalarize_dark_index + + +def test_scalarize_dark_index_accepts_length_one_numpy_arrays(): + assert _scalarize_dark_index(np.array([3])) == 3 + + +def test_scalarize_dark_index_accepts_numpy_scalars(): + assert _scalarize_dark_index(np.array(4)) == 4 + + +def test_scalarize_dark_index_accepts_uniform_length_many_arrays(): + assert _scalarize_dark_index(np.array([5, 5, 5])) == 5 + + +def test_scalarize_dark_index_rejects_nonuniform_arrays(): + with pytest.raises(ValueError, match="not uniform"): + _scalarize_dark_index(np.array([1, 2]))