|
| 1 | +# Copyright © UChicago Argonne LLC |
| 2 | +# See LICENSE file for details |
| 3 | +"""Synthetic result files: the datasets of a real multi-tau and two-time result, filled with a |
| 4 | +made-up scattering pattern and correlations, for tests, demos and continuous integration. |
| 5 | +
|
| 6 | +``pyxpcsviewer-synthetic FOLDER`` writes a small set; :func:`write_synthetic_result` writes one |
| 7 | +file. The files read back through :class:`~pyxpcsviewer.core.xpcs_file.XpcsFile` and every web |
| 8 | +route like a beamline result, only smaller. |
| 9 | +""" |
| 10 | + |
| 11 | +from __future__ import annotations |
| 12 | + |
| 13 | +import argparse |
| 14 | +import datetime as _dt |
| 15 | +import os |
| 16 | +import sys |
| 17 | +from pathlib import Path |
| 18 | + |
| 19 | +import h5py |
| 20 | +import numpy as np |
| 21 | + |
| 22 | +from .aps_8idi import key as hdf_key |
| 23 | + |
| 24 | +KEYS = hdf_key["nexus"] |
| 25 | + |
| 26 | + |
| 27 | +def _pattern(shape: tuple[int, int], center: tuple[float, float], rng: np.random.Generator) -> np.ndarray: |
| 28 | + """A SAXS-like intensity: a power-law decay from the beam center with a ring, a beam stop and noise.""" |
| 29 | + rows, cols = shape |
| 30 | + y = np.arange(rows)[:, None] - center[1] |
| 31 | + x = np.arange(cols)[None, :] - center[0] |
| 32 | + r = np.hypot(x, y) + 1e-3 |
| 33 | + intensity = 2.0e4 / (1 + (r / 12.0) ** 3.2) + 25.0 * np.exp(-((r - 0.42 * min(shape)) ** 2) / (2 * 9.0**2)) + 0.4 |
| 34 | + intensity *= 1 + 0.04 * rng.standard_normal(shape) |
| 35 | + intensity[r < 6] = 0.0 |
| 36 | + return intensity.astype(np.float32) |
| 37 | + |
| 38 | + |
| 39 | +def _g2_curves(t_el: np.ndarray, q: np.ndarray, rng: np.random.Generator) -> tuple[np.ndarray, np.ndarray]: |
| 40 | + """Stretched exponentials whose relaxation time falls with q: g2 (delays × bins) and its error.""" |
| 41 | + tau_q = 0.08 * (q[0] / q) |
| 42 | + g2 = 1.0 + 0.25 * np.exp(-2.0 * (t_el[:, None] / tau_q[None, :]) ** 0.9) |
| 43 | + err = 0.004 + 0.002 * (t_el[:, None] / t_el[-1]) |
| 44 | + g2 = g2 + rng.standard_normal(g2.shape) * err * 0.5 |
| 45 | + return g2.astype(np.float64), np.broadcast_to(err, g2.shape).astype(np.float64) |
| 46 | + |
| 47 | + |
| 48 | +def _c2_half(frames: int, tau_frames: float, rng: np.random.Generator) -> np.ndarray: |
| 49 | + """The upper triangle (with the diagonal) of a two-time correlation that decays with |t1 - t2|.""" |
| 50 | + i = np.arange(frames) |
| 51 | + lag = np.abs(i[:, None] - i[None, :]) |
| 52 | + c2 = 1.0 + 0.25 * np.exp(-2.0 * (lag / tau_frames) ** 0.9) + 0.01 * rng.standard_normal((frames, frames)) |
| 53 | + c2 = np.triu(c2) |
| 54 | + return c2.astype(np.float32) |
| 55 | + |
| 56 | + |
| 57 | +def write_synthetic_result( |
| 58 | + path: str | os.PathLike[str], |
| 59 | + *, |
| 60 | + multitau: bool = True, |
| 61 | + twotime: bool = True, |
| 62 | + shape: tuple[int, int] = (256, 240), |
| 63 | + frames: int = 400, |
| 64 | + nq: int = 8, |
| 65 | + nphi: int = 1, |
| 66 | + segments: int = 4, |
| 67 | + twotime_bins: int = 3, |
| 68 | + frame_time: float = 0.01, |
| 69 | + seed: int = 0, |
| 70 | +) -> Path: |
| 71 | + """Write one synthetic result file. |
| 72 | +
|
| 73 | + Args: |
| 74 | + path: Where to write (``.hdf``). |
| 75 | + multitau: Store a multi-tau analysis (g2, partial g2, delays). |
| 76 | + twotime: Store a two-time analysis (C2 maps of the first ``twotime_bins`` q-bins). |
| 77 | + shape: Detector rows, columns. |
| 78 | + frames: Frames of the acquisition (the C2 maps are ``frames`` square). |
| 79 | + nq: Dynamic q-bins (rings around the beam center). |
| 80 | + nphi: Static φ slices per q-bin. |
| 81 | + segments: Time segments of the partial g2 and the stability curves. |
| 82 | + twotime_bins: How many of the q-bins carry a C2 map. |
| 83 | + frame_time: Seconds per frame. |
| 84 | + seed: Of the noise. |
| 85 | + """ |
| 86 | + rng = np.random.default_rng(seed) |
| 87 | + rows, cols = shape |
| 88 | + center = (cols * 0.47, rows * 0.52) |
| 89 | + saxs_2d = _pattern(shape, center, rng) |
| 90 | + |
| 91 | + # ── the q-map: concentric rings as dynamic bins, the same rings cut in φ as static bins ── |
| 92 | + pixel_size = 75e-6 |
| 93 | + det_dist = 8.0 |
| 94 | + energy = 12.0 |
| 95 | + k0 = 2 * np.pi / (12.398 / energy) |
| 96 | + y = np.arange(rows)[:, None] - center[1] |
| 97 | + x = np.arange(cols)[None, :] - center[0] |
| 98 | + r = np.hypot(x, y) |
| 99 | + np.sin(np.arctan(r * pixel_size / det_dist)) * k0 |
| 100 | + phi = np.rad2deg(np.arctan2(y, x) * -1) |
| 101 | + r_min, r_max = 10.0, 0.45 * min(shape) |
| 102 | + edges = np.linspace(r_min, r_max, nq + 1) |
| 103 | + dqmap = np.zeros(shape, dtype=np.int32) |
| 104 | + for n in range(nq): |
| 105 | + dqmap[(r >= edges[n]) & (r < edges[n + 1])] = n + 1 |
| 106 | + phi_edges = np.linspace(-180, 180, nphi + 1) |
| 107 | + phi_index = np.clip(np.digitize(phi, phi_edges) - 1, 0, nphi - 1) |
| 108 | + sqmap = np.where(dqmap > 0, (dqmap - 1) * nphi + phi_index + 1, 0).astype(np.int32) |
| 109 | + mask = np.ones(shape, dtype=np.int32) |
| 110 | + mask[r < 6] = 0 |
| 111 | + mask[:, cols // 2 - 1 : cols // 2 + 1] = 0 # a module gap |
| 112 | + dqmap[mask == 0] = 0 |
| 113 | + sqmap[mask == 0] = 0 |
| 114 | + ring_r = 0.5 * (edges[:-1] + edges[1:]) |
| 115 | + dqlist = np.sin(np.arctan(ring_r * pixel_size / det_dist)) * k0 |
| 116 | + sqlist = dqlist.copy() |
| 117 | + dplist = np.array([0.0]) |
| 118 | + splist = 0.5 * (phi_edges[:-1] + phi_edges[1:]) if nphi > 1 else np.array([0.0]) |
| 119 | + dynamic_index_mapping = np.arange(nq) |
| 120 | + static_index_mapping = np.arange(nq * nphi) |
| 121 | + |
| 122 | + # ── 1-D: one value per static bin, the stability per segment ── |
| 123 | + counts = np.bincount(sqmap.ravel(), weights=saxs_2d.ravel().astype(np.float64), minlength=nq * nphi + 1)[1:] |
| 124 | + pixels = np.bincount(sqmap.ravel(), minlength=nq * nphi + 1)[1:] |
| 125 | + saxs_1d = (counts / np.maximum(pixels, 1)).astype(np.float64) |
| 126 | + drift = 1 + 0.03 * np.linspace(-1, 1, segments) |
| 127 | + iqp = (saxs_1d[None, :] * drift[:, None] * (1 + 0.01 * rng.standard_normal((segments, nq * nphi)))).astype( |
| 128 | + np.float64 |
| 129 | + ) |
| 130 | + |
| 131 | + # ── the intensity trace ── |
| 132 | + frame_index = np.arange(frames, dtype=np.float64) |
| 133 | + trace = float(saxs_2d[mask == 1].mean()) * ( |
| 134 | + 1 + 0.02 * rng.standard_normal(frames) + 0.01 * np.sin(frame_index / 25.0) |
| 135 | + ) |
| 136 | + int_t = np.stack([frame_index, trace]).astype(np.float64) |
| 137 | + |
| 138 | + # ── multi-tau ── |
| 139 | + delays = np.unique( |
| 140 | + np.concatenate([np.arange(1, 8), np.round(np.logspace(np.log10(8), np.log10(frames // 2), 24))]) |
| 141 | + ).astype(np.int64) |
| 142 | + t_el = delays * frame_time |
| 143 | + g2, g2_err = _g2_curves(t_el, dqlist, rng) |
| 144 | + g2_partial = np.stack([g2 * (1 + 0.004 * k * rng.standard_normal(g2.shape)) for k in range(segments)]) |
| 145 | + g2_partial_err = np.broadcast_to(g2_err, g2_partial.shape).copy() |
| 146 | + |
| 147 | + # ── two-time ── |
| 148 | + tau_frames = 0.08 / frame_time * (dqlist[0] / dqlist[:twotime_bins]) |
| 149 | + |
| 150 | + target = Path(path) |
| 151 | + target.parent.mkdir(parents=True, exist_ok=True) |
| 152 | + with h5py.File(target, "w") as f: |
| 153 | + |
| 154 | + def put(key: str, value: object, units: str | None = None) -> None: |
| 155 | + dataset = f.create_dataset(KEYS.get(key, key), data=value) |
| 156 | + if units: |
| 157 | + dataset.attrs["units"] = units |
| 158 | + |
| 159 | + put("mask", mask) |
| 160 | + # the q-map manager caches q-maps by this attribute: files of one geometry share one q-map |
| 161 | + f["/xpcs/qmap"].attrs["hash"] = f"synthetic-qmap-{rows}x{cols}-{nq}-{nphi}" |
| 162 | + put("dqmap", dqmap) |
| 163 | + put("sqmap", sqmap) |
| 164 | + put("dqlist", dqlist) |
| 165 | + put("sqlist", sqlist) |
| 166 | + put("dplist", dplist) |
| 167 | + put("splist", splist) |
| 168 | + put("bcx", float(center[0])) |
| 169 | + put("bcy", float(center[1])) |
| 170 | + put("det_dist", det_dist, "m") |
| 171 | + put("pixel_size", pixel_size, "m") |
| 172 | + put("static_index_mapping", static_index_mapping) |
| 173 | + put("dynamic_index_mapping", dynamic_index_mapping) |
| 174 | + put("static_num_pts", np.array([nq, nphi])) |
| 175 | + put("dynamic_num_pts", np.array([nq, 1])) |
| 176 | + put("map_names", np.array([b"q", b"phi"])) |
| 177 | + put("map_units", np.array(["Å⁻¹".encode(), b"deg"])) |
| 178 | + put("start_time", str(_dt.datetime(2026, 1, 1, 12, 0, 0)).encode()) |
| 179 | + put("saxs_2d", saxs_2d) |
| 180 | + put("saxs_1d", saxs_1d) |
| 181 | + put("Iqp", iqp) |
| 182 | + put("Int_t", int_t) |
| 183 | + put("t0", frame_time, "s") |
| 184 | + put("t1", frame_time * 0.9, "s") |
| 185 | + put("X_energy", energy, "keV") |
| 186 | + put("pix_dim_x", pixel_size, "m") |
| 187 | + put("pix_dim_y", pixel_size, "m") |
| 188 | + f.create_dataset("/entry/instrument/detector_1/description", data=b"synthetic detector") |
| 189 | + f.create_dataset("/entry/instrument/detector_1/x_pixels", data=cols) |
| 190 | + f.create_dataset("/entry/instrument/detector_1/y_pixels", data=rows) |
| 191 | + f.create_dataset("/entry/sample/name", data=b"synthetic sample") |
| 192 | + f.create_dataset("/entry/sample/temperature", data=295.0).attrs["units"] = "K" |
| 193 | + f.create_dataset("/entry/user/name", data=b"pyxpcsviewer-synthetic") |
| 194 | + if multitau: |
| 195 | + put("tau", delays) |
| 196 | + put("g2", g2) |
| 197 | + put("g2_err", g2_err) |
| 198 | + put("g2_partial", g2_partial) |
| 199 | + put("g2_partial_err", g2_partial_err) |
| 200 | + put("g2_partial_labels", np.array([f"{k + 1}".encode() for k in range(segments)])) |
| 201 | + put("stride_frame", 1) |
| 202 | + put("avg_frame", 1) |
| 203 | + f.create_dataset("/xpcs/multitau/config/analysis_type", data=b"Multitau") |
| 204 | + if twotime: |
| 205 | + prefix = KEYS["c2_prefix"] |
| 206 | + for n in range(twotime_bins): |
| 207 | + f.create_dataset(f"{prefix}/C2T_all_{n + 1:04d}", data=_c2_half(frames, float(tau_frames[n]), rng)) |
| 208 | + put("c2_processed_bins", np.arange(1, twotime_bins + 1)) |
| 209 | + lag = np.arange(frames, dtype=np.float64) |
| 210 | + c2_g2 = 1.0 + 0.25 * np.exp(-2.0 * (lag[:, None] / tau_frames[None, :]) ** 0.9) |
| 211 | + c2_g2[0] = 1.0 |
| 212 | + put("c2_g2", c2_g2.astype(np.float64)) |
| 213 | + seg_len = frames // segments |
| 214 | + seg_lag = np.arange(seg_len, dtype=np.float64) |
| 215 | + c2_g2_segments = 1.0 + 0.25 * np.exp(-2.0 * (seg_lag[:, None, None] / tau_frames[None, None, :]) ** 0.9) |
| 216 | + c2_g2_segments = c2_g2_segments * (1 + 0.005 * rng.standard_normal((seg_len, segments, twotime_bins))) |
| 217 | + put("c2_g2_segments", c2_g2_segments.astype(np.float64)) |
| 218 | + put("c2_delay", lag * frame_time) |
| 219 | + put("c2_stride_frame", 1) |
| 220 | + put("c2_avg_frame", 1) |
| 221 | + f.create_dataset("/xpcs/twotime/config/analysis_type", data=b"Twotime") |
| 222 | + return target |
| 223 | + |
| 224 | + |
| 225 | +def write_synthetic_folder(folder: str | os.PathLike[str], *, seed: int = 0) -> list[Path]: |
| 226 | + """A small folder of synthetic results: a multi-tau file, a two-time file and one with both.""" |
| 227 | + base = Path(folder) |
| 228 | + return [ |
| 229 | + write_synthetic_result(base / "S0001_synthetic_a0001_f000400_r00001_multitau.hdf", twotime=False, seed=seed), |
| 230 | + write_synthetic_result( |
| 231 | + base / "S0002_synthetic_a0002_f000400_r00001_twotime.hdf", multitau=False, seed=seed + 1 |
| 232 | + ), |
| 233 | + write_synthetic_result(base / "S0003_synthetic_a0003_f000400_r00001_both.hdf", seed=seed + 2), |
| 234 | + ] |
| 235 | + |
| 236 | + |
| 237 | +def main(argv: list[str] | None = None) -> int: |
| 238 | + parser = argparse.ArgumentParser( |
| 239 | + prog="pyxpcsviewer-synthetic", description="write synthetic XPCS result files for tests and demos" |
| 240 | + ) |
| 241 | + parser.add_argument("folder", help="where to write the files (created when missing)") |
| 242 | + parser.add_argument("--seed", type=int, default=0) |
| 243 | + args = parser.parse_args(argv) |
| 244 | + for path in write_synthetic_folder(args.folder, seed=args.seed): |
| 245 | + print(path, file=sys.stderr) |
| 246 | + return 0 |
| 247 | + |
| 248 | + |
| 249 | +if __name__ == "__main__": |
| 250 | + sys.exit(main()) |
0 commit comments