From 29a5fbd8e2f901f0d95ed68505416d17406edcc6 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 10 Jul 2026 08:51:51 +0000 Subject: [PATCH] Add petrolib.data_qc_io: cleaning, scaling, filtering, signal, synth, IO Final plan domain (LIBRARY_MERGE_PLAN.md data_qc_io). Package with six submodules, each grounded in a corpus survey of the srcYYYY_MM article forms and pre-verified bit-exact against the delegation targets: - clean: sentinels_to_nan, impute_gaps, outlier_mask (zscore/iqr/mad, sided), MAD despike (interp/median/nan), closure_residual, renormalize, relative_discrepancy (both corpus arrangements round identically) - scale: normalize_to_reference (Shier two-point, percentile or min/max endpoints), match_moments - filt: smooth (boxcar/gaussian/median), moving_stat (centered or causal; mean/std/var/mad/cv), tool_response (Gaussian vertical response), median_filter2d, window_features (edge-padded stacks), detect_bed_boundaries (|gradient|*sqrt(var) local maxima) - signal: snr_db (noise record or noise_std, VMD guard), stack_gain, block_stack, add_gaussian_noise (sigma = fraction * mean|x|) - synth: blocky_log, log_suite (verbatim CPA draw sequence), shifted_pair, gaussian_mixture_spectrum, sphere_pack_volume, disk_image (verbatim SEM draw sequence) - io: WellboreData (Bradley et al. 2023 universal JSON container, plus from_json) Deliberately not duplicated: unit conversions (petrolib.units), metrics and zscore/minmax (petrolib.ml_stats), xcorr/DTW alignment (petrolib.depth_matching), pay flags and net-to-gross (petrolib.porosity_lithology). 37/37 bit-exact pre-checks against article bodies; ruff + mypy clean; 252 tests pass (12 new). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01B2BX1XL8422gebCU1oVoQC --- petrolib/__init__.py | 1 + petrolib/data_qc_io/__init__.py | 40 +++++ petrolib/data_qc_io/clean.py | 174 +++++++++++++++++++++ petrolib/data_qc_io/filt.py | 203 ++++++++++++++++++++++++ petrolib/data_qc_io/io.py | 74 +++++++++ petrolib/data_qc_io/scale.py | 57 +++++++ petrolib/data_qc_io/signal.py | 86 +++++++++++ petrolib/data_qc_io/synth.py | 228 +++++++++++++++++++++++++++ tests/petrolib/test_data_qc_io.py | 246 ++++++++++++++++++++++++++++++ 9 files changed, 1109 insertions(+) create mode 100644 petrolib/data_qc_io/__init__.py create mode 100644 petrolib/data_qc_io/clean.py create mode 100644 petrolib/data_qc_io/filt.py create mode 100644 petrolib/data_qc_io/io.py create mode 100644 petrolib/data_qc_io/scale.py create mode 100644 petrolib/data_qc_io/signal.py create mode 100644 petrolib/data_qc_io/synth.py create mode 100644 tests/petrolib/test_data_qc_io.py diff --git a/petrolib/__init__.py b/petrolib/__init__.py index 9cf0d38..f0befd7 100644 --- a/petrolib/__init__.py +++ b/petrolib/__init__.py @@ -24,6 +24,7 @@ "borehole_image", "capillary_pressure", "constants", + "data_qc_io", "depth_correction", "depth_matching", "em_dielectric", diff --git a/petrolib/data_qc_io/__init__.py b/petrolib/data_qc_io/__init__.py new file mode 100644 index 0000000..5ceeac8 --- /dev/null +++ b/petrolib/data_qc_io/__init__.py @@ -0,0 +1,40 @@ +"""Log data QC, cleaning, filtering, synthetic data, and wellbore IO. + +The cross-cutting data-handling layer the articles re-implement around their +physics (see LIBRARY_MERGE_PLAN.md). Submodules: + +* :mod:`~petrolib.data_qc_io.clean` -- no-data sentinels, gap imputation, + outlier masks, despiking, and compositional closure checks. +* :mod:`~petrolib.data_qc_io.scale` -- reference-well (Shier) normalization + and moment matching. +* :mod:`~petrolib.data_qc_io.filt` -- smoothing, moving-window statistics, + Gaussian tool response, 2-D median filter, depth-window feature stacks, + and bed-boundary detection. +* :mod:`~petrolib.data_qc_io.signal` -- SNR, stacking gain, block stacking, + Gaussian noise injection. +* :mod:`~petrolib.data_qc_io.synth` -- synthetic blocky logs, log suites, + shifted curve pairs, Gaussian-mixture spectra, sphere packs, disk images. +* :mod:`~petrolib.data_qc_io.io` -- the Bradley et al. (2023) universal + wellbore-data JSON container. + +Related canon (deliberately NOT duplicated here): unit conversions live in +:mod:`petrolib.units`; error metrics and z-score / min-max feature scaling in +:mod:`petrolib.ml_stats`; depth alignment (cross-correlation, DTW, shifts) in +:mod:`petrolib.depth_matching`; pay flags, interval thickness, and +net-to-gross in :mod:`petrolib.porosity_lithology`. +""" + +from __future__ import annotations + +from . import clean, filt, io, scale, signal, synth +from .io import WellboreData + +__all__ = [ + "WellboreData", + "clean", + "filt", + "io", + "scale", + "signal", + "synth", +] diff --git a/petrolib/data_qc_io/clean.py b/petrolib/data_qc_io/clean.py new file mode 100644 index 0000000..40da113 --- /dev/null +++ b/petrolib/data_qc_io/clean.py @@ -0,0 +1,174 @@ +"""Log cleaning: sentinels, gap imputation, outliers, compositional closure. + +Conventions observed across the corpus and kept here: statistical spreads use +population std (``ddof=0``); comparison against outlier bounds is strict +(``>`` / ``<``); the robust MAD scale uses the Gaussian-consistency factor +1.4826. Sources: src2021_12/article01 (sentinels, z-score/IQR masks), +src2018_12/article11 (gap imputation, IQR), src2024_08/magnetic_permeability +(MAD despike), src2019_02/article4 (closure), src2023_08/article3 +(renormalize), src2016_02/article3 & src2017_10/article5 (discrepancy). +""" + +from __future__ import annotations + +from collections.abc import Sequence + +import numpy as np +from numpy.typing import ArrayLike, NDArray + +_Float = NDArray[np.float64] +_Bool = NDArray[np.bool_] + +#: Common LAS/DLIS no-data sentinel values. +SENTINELS = (-999.0, -999.25, -9999.0) + +#: Gaussian-consistency factor: sigma ~= 1.4826 * MAD for normal data. +MAD_TO_SIGMA = 1.4826 + + +def _arr(x: ArrayLike) -> _Float: + return np.asarray(x, np.float64) + + +def sentinels_to_nan( + x: ArrayLike, + sentinels: Sequence[float] = SENTINELS, + *, + rtol: float = 1e-5, + atol: float = 1e-8, +) -> _Float: + """Replace LAS/DLIS no-data sentinels (e.g. -999.25) with NaN. + + Matching uses ``np.isclose`` (defaults ``rtol=1e-5, atol=1e-8``) so float + round-trips of the sentinel still match. Source: src2021_12/article01. + """ + out = _arr(x).copy() + for s in sentinels: + out[np.isclose(out, s, rtol=rtol, atol=atol)] = np.nan + return out + + +def impute_gaps(x: ArrayLike, index: ArrayLike | None = None) -> _Float: + """Fill NaN runs by linear interpolation over the valid samples. + + ``index`` gives the sample coordinates (default ``arange``); leading and + trailing NaNs are held flat at the nearest valid value (``np.interp`` + endpoint behavior). Source: src2018_12/article11. + """ + out = _arr(x).copy() + idx = np.arange(out.size) if index is None else _arr(index) + good = ~np.isnan(out) + if not good.any(): + raise ValueError("impute_gaps: input has no finite values to interpolate from") + out[~good] = np.interp(idx[~good], idx[good], out[good]) + return out + + +def outlier_mask( + x: ArrayLike, + method: str = "zscore", + *, + threshold: float = 3.0, + k: float = 1.5, + side: str = "both", +) -> _Bool: + """Boolean outlier mask by z-score, IQR whiskers, or robust MAD. + + ``'zscore'``: ``|z| > threshold`` with population std (guarded ``>1e-12``). + ``'iqr'``: outside ``(Q1 - k*IQR, Q3 + k*IQR)``. + ``'mad'``: further than ``threshold * MAD * 1.4826`` from the median. + ``side`` restricts to ``'upper'``/``'lower'`` exceedances (the + ferromagnetic-contaminant screen is upper-only). Sources: + src2021_12/article01, src2018_12/article11, src2024_08/magnetic_permeability. + """ + if side not in {"both", "upper", "lower"}: + raise ValueError(f"outlier_mask: unknown side {side!r}") + arr = _arr(x) + if method == "zscore": + s = arr.std() + z = (arr - arr.mean()) / (s if s > 1e-12 else 1.0) + if side == "both": + return np.asarray(np.abs(z) > threshold) + return np.asarray(z > threshold if side == "upper" else z < -threshold) + if method == "iqr": + q1, q3 = np.percentile(arr, [25, 75]) + iqr = q3 - q1 + lo, hi = q1 - k * iqr, q3 + k * iqr + if side == "both": + return np.asarray((arr < lo) | (arr > hi)) + return np.asarray(arr > hi if side == "upper" else arr < lo) + if method == "mad": + med = np.median(arr) + h = threshold * np.median(np.abs(arr - med)) * MAD_TO_SIGMA + if side == "both": + return np.asarray((arr < med - h) | (arr > med + h)) + return np.asarray(arr > med + h if side == "upper" else arr < med - h) + raise ValueError(f"outlier_mask: unknown method {method!r}") + + +def despike( + x: ArrayLike, + *, + threshold: float = 3.0, + side: str = "both", + replace: str = "interp", +) -> _Float: + """Remove spikes flagged by the robust MAD criterion. + + Spikes are samples further than ``threshold * MAD * 1.4826`` from the + median (see :func:`outlier_mask`). ``replace='interp'`` fills them by + linear interpolation from the surviving samples, ``'median'`` with the + global median (the ferromagnetic-contaminant convention), ``'nan'`` with + NaN. Source: src2024_08/magnetic_permeability (one-sided, median). + """ + out = _arr(x).copy() + spikes = outlier_mask(out, "mad", threshold=threshold, side=side) + if not spikes.any(): + return out + if replace == "interp": + idx = np.arange(out.size) + good = ~spikes + out[spikes] = np.interp(idx[spikes], idx[good], out[good]) + elif replace == "median": + out[spikes] = np.median(_arr(x)) + elif replace == "nan": + out[spikes] = np.nan + else: + raise ValueError(f"despike: unknown replace {replace!r}") + return out + + +def closure_residual(*fracs: ArrayLike, target: float = 1.0) -> _Float: + """Closure residual ``sum(fracs) - target`` (e.g. ``Sw+So+Sg-1``). + + Zero means the fractions close exactly. Source: src2019_02/article4 + (saturation closure); the oxide-closure factor lives in + ``petrolib.geochem_fluids``. + """ + if not fracs: + raise ValueError("closure_residual: at least one fraction is required") + total = _arr(fracs[0]) + for f in fracs[1:]: + total = total + _arr(f) + return np.asarray(total - target) + + +def renormalize(fracs: ArrayLike, axis: int = -1) -> _Float: + """Rescale fractions to sum to one along ``axis`` (compositional closure). + + A zero sum yields NaN/inf rather than a guessed composition. Sources: + src2023_08/article3 (mineral volumes), src2014_10/article2 (oxides, via + ``petrolib.geochem_fluids.core_geochem.oxide_closure``). + """ + f = _arr(fracs) + return np.asarray(f / np.sum(f, axis=axis, keepdims=True)) + + +def relative_discrepancy(a: ArrayLike, b: ArrayLike) -> _Float: + """Relative discrepancy ``2|a-b| / (a+b)`` (difference over the mean). + + Both corpus arrangements (``2|a-b|/(a+b)`` and ``|a-b|/(0.5(a+b))``) + round identically in IEEE arithmetic. Sources: src2016_02/article3 + (deviation indicator, Eq. 6), src2017_10/article5. + """ + return np.asarray(2.0 * np.abs(_arr(a) - _arr(b)) / (_arr(a) + _arr(b))) diff --git a/petrolib/data_qc_io/filt.py b/petrolib/data_qc_io/filt.py new file mode 100644 index 0000000..8d3dd6f --- /dev/null +++ b/petrolib/data_qc_io/filt.py @@ -0,0 +1,203 @@ +"""Depth-series filtering: smoothing, moving statistics, tool response, +2-D median filter, window feature stacks, bed-boundary detection. + +Edge conventions differ per idiom and are preserved: kernel smoothing +(boxcar/Gaussian) zero-pads via ``np.convolve(mode='same')``; the median +smoother and the feature stack edge-pad (``np.pad(mode='edge')``, i.e. +scipy's ``'nearest'``); moving statistics shrink the window at the edges. +All spreads are population (``ddof=0``). +""" + +from __future__ import annotations + +import numpy as np +from numpy.typing import ArrayLike, NDArray + +_Float = NDArray[np.float64] + +_MOVING_STATS = ("mean", "std", "var", "mad", "cv") + + +def _arr(x: ArrayLike) -> _Float: + return np.asarray(x, np.float64) + + +def smooth( + x: ArrayLike, + window: int = 5, + kind: str = "boxcar", + *, + sigma: float | None = None, +) -> _Float: + """Smooth a 1-D curve with a boxcar, Gaussian, or median filter. + + ``'boxcar'``: centered moving average, ``np.convolve(mode='same')`` + (zero-padded edges); ``window<=1`` returns a copy. ``'gaussian'``: + normalized Gaussian kernel of standard deviation ``sigma`` samples + (kernel length ``max(3, round(6*sigma)|1)``), same convolution. + ``'median'``: running median over an odd ``window`` with edge-replicated + padding. Sources: src2023_10/article_03 (boxcar), src2026_04/a04 + (moving_average), src2024_06/article2 (median, mode='nearest'). + """ + arr = _arr(x) + if kind == "boxcar": + if window <= 1: + return arr.copy() + k = np.ones(window) / window + return np.asarray(np.convolve(arr, k, mode="same")) + if kind == "gaussian": + if sigma is None: + raise ValueError("smooth: sigma is required for kind='gaussian'") + n = max(3, int(round(6.0 * sigma)) | 1) + t = np.arange(n) - n // 2 + g = np.exp(-0.5 * (t / sigma) ** 2) + g /= g.sum() + return np.asarray(np.convolve(arr, g, mode="same")) + if kind == "median": + if window < 1 or window % 2 == 0: + raise ValueError("smooth: kind='median' requires an odd window") + pad = window // 2 + padded = np.pad(arr, pad, mode="edge") + windows = np.lib.stride_tricks.sliding_window_view(padded, window) + return np.asarray(np.median(windows, axis=-1)) + raise ValueError(f"smooth: unknown kind {kind!r}") + + +def moving_stat( + x: ArrayLike, + window: int, + stat: str = "mean", + *, + center: bool = True, +) -> _Float: + """Moving-window statistic with shrinking edge windows. + + ``center=True`` uses the symmetric window ``x[i-window//2 : i+window//2+1]`` + (clamped at the edges); ``center=False`` uses the causal/trailing window + ``x[i-window+1 : i+1]``. ``stat`` is ``'mean'``, ``'std'``, ``'var'`` + (all ``ddof=0``), ``'mad'`` (mean absolute deviation about the window + mean), or ``'cv'`` (``std/mean``, 0 where the mean is 0). Sources: + src2020_10/article3 (rolling_cv, half-width convention), src2024_10/ + ml_permeability (centered mean/var), src2026_04/a04 (causal std/msd). + """ + if stat not in _MOVING_STATS: + raise ValueError(f"moving_stat: unknown stat {stat!r}") + arr = _arr(x) + n = arr.size + half = window // 2 + out = np.zeros(n) + for i in range(n): + if center: + seg = arr[max(0, i - half) : min(n, i + half + 1)] + else: + seg = arr[max(0, i - window + 1) : i + 1] + if stat == "mean": + out[i] = np.mean(seg) + elif stat == "std": + out[i] = np.std(seg) + elif stat == "var": + out[i] = np.var(seg) + elif stat == "mad": + out[i] = np.mean(np.abs(seg - np.mean(seg))) + else: # cv + m = seg.mean() + out[i] = seg.std() / m if m != 0 else 0.0 + return out + + +def tool_response(x: ArrayLike, dz: float, fwhm: float) -> _Float: + """Convolve a fine-scale curve with a Gaussian tool vertical response. + + ``sigma = fwhm / (2*sqrt(2*ln 2))``; the kernel spans ~6 sigma with an + odd length (``max(3, round(6*sigma/dz)|1)``), is unit-normalized, and is + applied with ``np.convolve(mode='same')``. This is the log-simulation / + upscaling operator for a tool of aperture ``fwhm`` on a grid of spacing + ``dz`` (same length units). Source: src2023_10/article_09 + (vertical_response). + """ + arr = _arr(x) + sigma = fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0))) + n = max(3, int(round(6 * sigma / dz)) | 1) + t = (np.arange(n) - n // 2) * dz + g = np.exp(-0.5 * (t / sigma) ** 2) + g /= g.sum() + return np.asarray(np.convolve(arr, g, mode="same")) + + +def median_filter2d(img: ArrayLike, size: int = 3) -> _Float: + """Pure-numpy 2-D median filter with edge-replicated padding. + + Builds the ``size*size`` stack of shifted views of the padded image and + takes the median across it. Source: src2023_06/article9 + (median_filter_3x3), generalized to odd ``size``. + """ + arr = _arr(img) + if arr.ndim != 2: + raise ValueError("median_filter2d: img must be 2-D") + if size < 1 or size % 2 == 0: + raise ValueError("median_filter2d: size must be odd") + pad = size // 2 + padded = np.pad(arr, pad, mode="edge") + h, w = arr.shape + stack = np.empty((size * size, h, w)) + k = 0 + for dy in range(size): + for dx in range(size): + stack[k] = padded[dy : dy + h, dx : dx + w] + k += 1 + return np.asarray(np.median(stack, axis=0)) + + +def window_features(x: ArrayLike, window: int = 11) -> _Float: + """Edge-padded sliding-window feature stack for depth-indexed ML. + + A 1-D curve of length ``n`` becomes ``(n, window)`` — row ``i`` is the + window centered on sample ``i``, with ``np.pad(mode='edge')`` handling + the ends so every depth keeps a full window. A 2-D ``(n, d)`` matrix + stacks each column in turn to ``(n, d*window)``. ``window`` must be odd. + Sources: src2023_04/article10 & article11 (induction deconvolution + feature stacks). + """ + if window < 1 or window % 2 == 0: + raise ValueError("window_features: window must be odd") + arr = _arr(x) + if arr.ndim == 1: + pad = window // 2 + padded = np.pad(arr, pad, mode="edge") + return np.lib.stride_tricks.sliding_window_view(padded, window).copy() + if arr.ndim == 2: + return np.asarray( + np.hstack([window_features(arr[:, j], window) for j in range(arr.shape[1])]) + ) + raise ValueError("window_features: x must be 1-D or 2-D") + + +def detect_bed_boundaries( + curve: ArrayLike, + window: int = 11, + threshold: float | None = None, +) -> NDArray[np.int_]: + """Bed-boundary indices from a derivative-times-local-variability metric. + + ``metric = |gradient(curve)| * sqrt(moving var)`` (variance over the + centered ``window``, edges left 0); boundaries are the local maxima of + the metric above ``threshold`` (default adaptive ``mean + 2*std``). + Source: src2023_08/article5 (detect_bed_boundaries); simpler + gradient-threshold and lag-variance variants exist in src2024_10/ + thin_bed_nmr and src2024_06/article2. + """ + arr = _arr(curve) + n = arr.size + half = window // 2 + deriv = np.gradient(arr) + variance = np.zeros(n) + for i in range(half, n - half): + variance[i] = np.var(arr[i - half : i + half + 1]) + metric = np.abs(deriv) * np.sqrt(variance) + if threshold is None: + threshold = float(np.mean(metric) + 2.0 * np.std(metric)) + boundaries = [] + for i in range(1, n - 1): + if metric[i] > threshold and metric[i] >= metric[i - 1] and metric[i] >= metric[i + 1]: + boundaries.append(i) + return np.array(boundaries, dtype=int) diff --git a/petrolib/data_qc_io/io.py b/petrolib/data_qc_io/io.py new file mode 100644 index 0000000..3e1d923 --- /dev/null +++ b/petrolib/data_qc_io/io.py @@ -0,0 +1,74 @@ +"""Universal wellbore-data container (Bradley et al. 2023). + +A minimal JSON-based format for storing and exchanging wellbore data: +hierarchical metadata plus named n-dimensional channel arrays with units and +axis labels, covering simple 1-D logs through multi-dimensional measurements +such as ultradeep azimuthal resistivity. Source: src2023_12/ +bradley_wellbore_format.py (Petrophysics 64(6):823-836). +""" + +from __future__ import annotations + +import json +from typing import Any + +import numpy as np +from numpy.typing import ArrayLike + + +class WellboreData: + """Hierarchical metadata + n-D channels with units and axes. + + ``meta`` holds ``well_name``, ``uwi``, and the format ``version``; + ``channels`` maps name -> ``{"data", "units", "axes", "shape"}``. + Round-trip via ``to_json``/``from_json`` (or ``to_dict``/``from_dict``); + ``shape`` is recomputed on reload. + """ + + def __init__(self, well_name: str, uwi: str | None = None) -> None: + self.meta: dict[str, Any] = {"well_name": well_name, "uwi": uwi, "version": "1.0"} + self.channels: dict[str, dict[str, Any]] = {} + + def add_channel(self, name: str, data: ArrayLike, units: str, axes: list[str]) -> None: + """Store an n-D channel; ``axes`` must name every array dimension.""" + arr = np.asarray(data) + if len(axes) != arr.ndim: + raise ValueError(f"axes count {len(axes)} != ndim {arr.ndim}") + self.channels[name] = { + "data": arr, + "units": units, + "axes": axes, + "shape": arr.shape, + } + + def to_dict(self) -> dict[str, Any]: + """JSON-serializable dict; channel arrays become nested lists.""" + return { + "meta": self.meta, + "channels": { + n: {**{k: v for k, v in c.items() if k != "data"}, "data": c["data"].tolist()} + for n, c in self.channels.items() + }, + } + + def to_json(self, path: str | None = None) -> str: + """Serialize to a JSON string (``indent=2``); optionally write ``path``.""" + s = json.dumps(self.to_dict(), indent=2) + if path: + with open(path, "w") as fh: + fh.write(s) + return s + + @classmethod + def from_dict(cls, d: dict[str, Any]) -> WellboreData: + """Rebuild from :meth:`to_dict` output (or parsed JSON).""" + w = cls(d["meta"]["well_name"], d["meta"].get("uwi")) + w.meta = d["meta"] + for n, c in d["channels"].items(): + w.add_channel(n, np.array(c["data"]), c["units"], c["axes"]) + return w + + @classmethod + def from_json(cls, s: str) -> WellboreData: + """Rebuild from a JSON string produced by :meth:`to_json`.""" + return cls.from_dict(json.loads(s)) diff --git a/petrolib/data_qc_io/scale.py b/petrolib/data_qc_io/scale.py new file mode 100644 index 0000000..fc0acb0 --- /dev/null +++ b/petrolib/data_qc_io/scale.py @@ -0,0 +1,57 @@ +"""Curve normalization against a reference well or target moments. + +Only the reference-based normalizations live here; plain feature scaling is +already canon in :mod:`petrolib.ml_stats` (``zscore``, ``minmax``, +``affine_rescale``). Population std (``ddof=0``) throughout, matching every +article copy. +""" + +from __future__ import annotations + +import numpy as np +from numpy.typing import ArrayLike, NDArray + +_Float = NDArray[np.float64] + + +def _arr(x: ArrayLike) -> _Float: + return np.asarray(x, np.float64) + + +def normalize_to_reference( + x: ArrayLike, + ref_lo: float, + ref_hi: float, + *, + in_lo: float | None = None, + in_hi: float | None = None, + pct: tuple[float, float] | None = (5.0, 95.0), +) -> _Float: + """Two-point (Shier 2004) normalization onto reference-well endpoints. + + ``out = ref_lo + (ref_hi - ref_lo) * (x - in_lo) / (in_hi - in_lo)``. + + ``in_lo``/``in_hi`` are the input-well endpoints; when omitted they are + taken as the ``pct`` percentiles of ``x`` (the Shier 5th/95th convention), + or the min/max when ``pct=None`` (the image-histogram convention). No + clipping is applied. Sources: src2021_12/article01 (normalize_reference), + src2016_12/article6 (histogram_normalize), src2015_04/article4 + (histogram_scale). + """ + v = _arr(x) + if in_lo is None: + in_lo = float(np.percentile(v, pct[0])) if pct is not None else float(np.min(v)) + if in_hi is None: + in_hi = float(np.percentile(v, pct[1])) if pct is not None else float(np.max(v)) + return np.asarray(ref_lo + (ref_hi - ref_lo) * (v - in_lo) / (in_hi - in_lo)) + + +def match_moments(x: ArrayLike, target_mean: float, target_std: float) -> _Float: + """Affine transform of ``x`` to the target mean and (population) std. + + ``out = (x - mean(x)) / std(x) * target_std + target_mean``; a constant + input divides by zero rather than guessing a scale. Source: + src2016_12/article6 (affine_normalize). + """ + v = _arr(x) + return np.asarray((v - v.mean()) / v.std() * target_std + target_mean) diff --git a/petrolib/data_qc_io/signal.py b/petrolib/data_qc_io/signal.py new file mode 100644 index 0000000..3d01dab --- /dev/null +++ b/petrolib/data_qc_io/signal.py @@ -0,0 +1,86 @@ +"""Signal-to-noise, stacking, and controlled noise injection. + +SNR is the power ratio ``10*log10(P_signal/P_noise)``; stacking of ``n`` +repeats improves amplitude SNR by ``sqrt(n)``. +""" + +from __future__ import annotations + +import numpy as np +from numpy.typing import ArrayLike, NDArray + +_Float = NDArray[np.float64] + + +def _arr(x: ArrayLike) -> _Float: + return np.asarray(x, np.float64) + + +def snr_db( + signal: ArrayLike, + *, + noise: ArrayLike | None = None, + noise_std: float | None = None, + guard: float = 0.0, +) -> float: + """Signal-to-noise ratio in dB: ``10*log10(mean(s^2)/P_noise)``. + + Give either a noise record/residual (``noise``, power ``mean(n^2)+guard``) + or a known standard deviation (``noise_std``, power ``noise_std**2``). + Zero noise power returns ``inf``; ``guard`` reproduces the additive + ``+1e-30`` convention of the VMD paper. Sources: src2025_08/ + fiber_optics_sensing (compute_snr), src2025_08/anomaly_detection_vmd. + """ + sig_power = np.mean(_arr(signal) ** 2) + if noise_std is not None and noise is None: + noise_power = float(noise_std) ** 2 + elif noise is not None and noise_std is None: + noise_power = float(np.mean(_arr(noise) ** 2)) + guard + else: + raise ValueError("snr_db: give exactly one of noise or noise_std") + if noise_power == 0: + return float("inf") + return float(10.0 * np.log10(sig_power / noise_power)) + + +def stack_gain(n: ArrayLike) -> _Float: + """Amplitude SNR gain ``sqrt(n)`` from stacking ``n`` repeat measurements. + + Sources: src2021_04/article5 (stacking_snr), src2025_12/carbon13_mr. + """ + return np.asarray(np.sqrt(_arr(n))) + + +def block_stack(x: ArrayLike, n: int, axis: int = -1) -> _Float: + """Block-mean stacking: average consecutive groups of ``n`` along ``axis``. + + The output length is ``size // n``; a trailing partial block is dropped. + Source: src2025_08/fiber_optics_sensing (temporal_stack). + """ + if n < 1: + raise ValueError("block_stack: n must be >= 1") + arr = _arr(x) + moved = np.moveaxis(arr, axis, -1) + n_out = moved.shape[-1] // n + trimmed = moved[..., : n_out * n] + stacked = trimmed.reshape(*moved.shape[:-1], n_out, n).mean(axis=-1) + return np.asarray(np.moveaxis(stacked, -1, axis)) + + +def add_gaussian_noise( + x: ArrayLike, + sigma_fraction: float, + rng: np.random.Generator | int | None = None, +) -> _Float: + """Add zero-mean Gaussian noise with ``sigma = fraction * mean(|x|)``. + + The scale is a single scalar tied to the curve's mean absolute level. + ``rng=None`` uses the deterministic ``default_rng(0)`` of the source + article; pass a seed or Generator to vary. Source: src2021_12/article01 + (Eqs. 9-10). + """ + arr = _arr(x) + if not isinstance(rng, np.random.Generator): + rng = np.random.default_rng(0 if rng is None else rng) + sigma = sigma_fraction * np.abs(arr).mean() + return np.asarray(arr + rng.normal(0.0, sigma, size=arr.shape)) diff --git a/petrolib/data_qc_io/synth.py b/petrolib/data_qc_io/synth.py new file mode 100644 index 0000000..33d993f --- /dev/null +++ b/petrolib/data_qc_io/synth.py @@ -0,0 +1,228 @@ +"""Synthetic test data: blocky logs, log suites, shifted pairs, spectra, +sphere packs, and disk images. + +Generators are deterministic by default (``rng=None`` seeds ``default_rng`` +with the source article's seed) so examples and tests reproduce; pass a seed +or ``np.random.Generator`` to vary. +""" + +from __future__ import annotations + +import numpy as np +from numpy.typing import ArrayLike, NDArray + +_Float = NDArray[np.float64] +_Int = NDArray[np.int_] + +#: Curve order of :func:`log_suite` columns. +LOG_SUITE_CURVES = ("GR", "RD", "RHOB", "DT", "NPHI") + +#: Per-facies (rows) mean levels for the :data:`LOG_SUITE_CURVES` (columns): +#: granitic, shale, fractured granite, weathered. Source: src2022_08/article5. +LOG_SUITE_FACIES_PROPS = np.array( + [ + [22.0, 80.0, 2.45, 75.0, 0.08], + [75.0, 10.0, 2.20, 95.0, 0.30], + [35.0, 200.0, 2.60, 65.0, 0.05], + [50.0, 30.0, 2.40, 80.0, 0.15], + ] +) + + +def _arr(x: ArrayLike) -> _Float: + return np.asarray(x, np.float64) + + +def _rng(rng: np.random.Generator | int | None, default_seed: int = 0) -> np.random.Generator: + if isinstance(rng, np.random.Generator): + return rng + return np.random.default_rng(default_seed if rng is None else rng) + + +def blocky_log( + n: int, + n_beds: int = 8, + *, + level_range: tuple[float, float] = (20.0, 150.0), + noise: float = 2.0, + rng: np.random.Generator | int | None = None, +) -> tuple[_Float, _Int]: + """Blocky (square) log of ``n_beds`` constant levels plus Gaussian noise. + + Returns ``(curve, bounds)`` where ``bounds`` are the ``n_beds - 1`` + interior boundary indices, drawn without replacement away from the ends + (margin ``max(1, n//15)``). Levels are uniform on ``level_range``. + Composed from the corpus generators: random breakpoints per + src2022_08/article5, curve+bounds return per src2024_06/article2, + ``noise=2.0`` per src2022_02/article2. + """ + if n_beds < 1: + raise ValueError("blocky_log: n_beds must be >= 1") + gen = _rng(rng) + margin = max(1, n // 15) + interior = np.arange(margin, n - margin) + if interior.size < n_beds - 1: + raise ValueError("blocky_log: n too small for n_beds") + bounds = np.sort(gen.choice(interior, n_beds - 1, replace=False)).astype(int) + edges = np.concatenate(([0], bounds, [n])) + levels = gen.uniform(level_range[0], level_range[1], n_beds) + curve = np.zeros(n) + for a, b, mu in zip(edges[:-1], edges[1:], levels, strict=True): + curve[a:b] = mu + return curve + gen.normal(0.0, noise, n), bounds + + +def log_suite( + n: int = 600, + n_facies: int = 4, + *, + rng: np.random.Generator | int | None = None, +) -> tuple[dict[str, _Float], _Int, _Int]: + """Five-log synthetic well (GR, RD, RHOB, DT, NPHI) with facies blocks. + + Returns ``(logs, boundaries, facies_labels)``: a dict of the + :data:`LOG_SUITE_CURVES`, the interior breakpoint indices, and the + per-sample facies index. Each block draws each curve from + ``N(mu, 0.05*|mu| + 1e-3)`` with ``mu`` from + :data:`LOG_SUITE_FACIES_PROPS`; breakpoints are drawn from + ``arange(40, n-40)`` without replacement. Source: src2022_08/article5 + (synth_log_suite, verbatim draw sequence). + """ + if not 1 <= n_facies <= len(LOG_SUITE_FACIES_PROPS): + raise ValueError(f"log_suite: n_facies must be in 1..{len(LOG_SUITE_FACIES_PROPS)}") + gen = _rng(rng) + bp_true = sorted(gen.choice(np.arange(40, n - 40), n_facies - 1, replace=False)) + bp_full = [0, *bp_true, n] + facies_labels = np.zeros(n, dtype=int) + out = np.zeros((n, 5)) + for k, (a, b) in enumerate(zip(bp_full[:-1], bp_full[1:], strict=True)): + fac = k % n_facies + facies_labels[a:b] = fac + for j in range(5): + mu = LOG_SUITE_FACIES_PROPS[fac, j] + sigma = 0.05 * abs(mu) + 1e-3 + out[a:b, j] = gen.normal(mu, sigma, b - a) + logs = {name: out[:, j] for j, name in enumerate(LOG_SUITE_CURVES)} + return logs, np.array(bp_true, dtype=int), facies_labels + + +def shifted_pair( + n: int, + shift: int, + *, + noise: float = 0.0, + rng: np.random.Generator | int | None = None, +) -> tuple[_Float, _Float]: + """Reference curve and a copy lagged by ``shift`` samples. + + Both come from one multi-harmonic base with ``other[i] = ref[i + shift]`` + (noise-free) — aligning ``other`` to ``ref`` (e.g. with + ``petrolib.depth_matching.xcorr_shift``) recovers ``shift``. ``noise`` + adds independent Gaussian noise to each curve. Pattern: + src2024_10/rddtw_depth_matching (integer-shift core/log pair), + src2022_02/article2 (rigid-lag windows). + """ + s = int(shift) + m = n + abs(s) + z = np.linspace(0.0, 4.0 * np.pi, m) + base = np.sin(z) + 0.5 * np.sin(3.0 * z) + 0.3 * np.cos(11.0 * z) + ref = base[max(-s, 0) : max(-s, 0) + n] + other = base[max(s, 0) : max(s, 0) + n] + if noise > 0.0: + gen = _rng(rng) + ref = ref + noise * gen.standard_normal(n) + other = other + noise * gen.standard_normal(n) + return np.asarray(ref), np.asarray(other) + + +def gaussian_mixture_spectrum( + axis: ArrayLike, + centres: ArrayLike, + amps: ArrayLike, + widths: ArrayLike, + *, + noise: float = 0.0, + log_axis: bool = True, + rng: np.random.Generator | int | None = None, +) -> _Float: + """Sum-of-Gaussians spectrum (synthetic NMR T2 / pore-size distributions). + + Each component contributes ``amp * exp(-0.5*((u - c)/width)**2)`` where + ``u`` is ``log10(axis)`` and ``c`` is ``log10(centre)`` when + ``log_axis=True`` (widths in decades), or the raw values otherwise — + the standard-deviation form of src2026_06/a03 (gaussian_t2_model) and + src2022_06/article1. ``noise>0`` adds Gaussian noise and clips at zero + (the src2022_06/article9 convention). + """ + ax = _arr(axis) + u = np.log10(ax) if log_axis else ax + out = np.zeros_like(u) + comps = zip(_arr(centres).ravel(), _arr(amps).ravel(), _arr(widths).ravel(), strict=True) + for c, a, w in comps: + cc = np.log10(c) if log_axis else c + out += a * np.exp(-0.5 * ((u - cc) / w) ** 2) + if noise > 0.0: + out = out + noise * _rng(rng).standard_normal(u.size) + out = np.clip(out, 0.0, None) + return out + + +def sphere_pack_volume( + shape: tuple[int, int, int], + r_range: tuple[int, int], + target_phi: float, + *, + rng: np.random.Generator | int | None = None, + max_grains: int = 100_000, +) -> NDArray[np.uint8]: + """Random sphere pack on a voxel grid (0 = pore, 1 = grain). + + Spheres with integer voxel radii uniform on ``r_range`` (inclusive) are + stamped at uniform random centers until the pore fraction drops to + ``target_phi`` (or ``max_grains`` is hit). Grain fill uses the + squared-distance broadcast of src2022_12/article7 (make_sand_pack). + """ + gen = _rng(rng) + nx, ny, nz = shape + cube = np.zeros(shape, dtype=np.uint8) + x = np.arange(nx)[:, None, None] + y = np.arange(ny)[None, :, None] + z = np.arange(nz)[None, None, :] + r_lo, r_hi = r_range + for _ in range(max_grains): + if 1.0 - cube.mean() <= target_phi: + break + cx = int(gen.integers(0, nx)) + cy = int(gen.integers(0, ny)) + cz = int(gen.integers(0, nz)) + r = int(gen.integers(r_lo, r_hi + 1)) + d2 = (x - cx) ** 2 + (y - cy) ** 2 + (z - cz) ** 2 + cube[d2 <= r * r] = 1 + return cube + + +def disk_image( + size: int, + n_features: int, + *, + rng: np.random.Generator | int | None = None, +) -> tuple[_Float, NDArray[np.bool_]]: + """Grey image with darker random disks; returns ``(image, mask)``. + + Background ``0.6 + 0.05*N(0,1)``, ``n_features`` disks of radius 2-7 px + at uniform centers, disk pixels ``0.15 + 0.05*N(0,1)``, clipped to + [0, 1] — a synthetic SEM/pore image. Source: src2024_04/ + chen_sem_pore_segmentation (synthetic_sem_image, verbatim draw + sequence; that article seeds ``default_rng(4)``). + """ + gen = _rng(rng) + img = 0.6 + 0.05 * gen.standard_normal((size, size)) + mask = np.zeros((size, size), dtype=bool) + yy, xx = np.mgrid[:size, :size] + for _ in range(n_features): + cy, cx = gen.integers(0, size, size=2) + r = gen.integers(2, 8) + m = (yy - cy) ** 2 + (xx - cx) ** 2 <= r * r + mask |= m + img[mask] = 0.15 + 0.05 * gen.standard_normal(mask.sum()) + return np.clip(img, 0, 1), mask diff --git a/tests/petrolib/test_data_qc_io.py b/tests/petrolib/test_data_qc_io.py new file mode 100644 index 0000000..219256f --- /dev/null +++ b/tests/petrolib/test_data_qc_io.py @@ -0,0 +1,246 @@ +"""Tests for petrolib.data_qc_io: cleaning, scaling, filtering, signal, +synthetic data, and the wellbore JSON container -- known values, article +anchors, and error paths.""" + +from __future__ import annotations + +import sys +from pathlib import Path + +import numpy as np +import pytest + +REPO_ROOT = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(REPO_ROOT)) + +from petrolib import data_qc_io as dq # noqa: E402 + +# --- clean ---------------------------------------------------------------------- + + +def test_sentinels_and_impute() -> None: + x = np.array([1.0, -999.25, 3.0, -9999.0, 5.0]) + out = dq.clean.sentinels_to_nan(x) + assert np.isnan(out[[1, 3]]).all() and np.array_equal(out[[0, 2, 4]], [1.0, 3.0, 5.0]) + # linear interpolation over the gap, flat extrapolation at the ends + filled = dq.clean.impute_gaps(np.array([np.nan, 1.0, np.nan, 3.0, np.nan])) + np.testing.assert_allclose(filled, [1.0, 1.0, 2.0, 3.0, 3.0]) + # explicit non-uniform index weights the interpolation + filled = dq.clean.impute_gaps(np.array([0.0, np.nan, 4.0]), index=[0.0, 3.0, 4.0]) + np.testing.assert_allclose(filled, [0.0, 3.0, 4.0]) + with pytest.raises(ValueError, match="no finite"): + dq.clean.impute_gaps(np.array([np.nan, np.nan])) + + +def test_outlier_mask_and_despike() -> None: + x = np.concatenate([np.zeros(50), [100.0]]) + for method in ("zscore", "iqr", "mad"): + mask = dq.clean.outlier_mask(x, method) + assert mask[-1] and mask[:50].sum() == 0 + # one-sided: an upper screen never flags low values + lo = np.concatenate([np.ones(50), [-100.0]]) + assert not dq.clean.outlier_mask(lo, "mad", side="upper").any() + assert dq.clean.outlier_mask(lo, "mad", side="lower")[-1] + with pytest.raises(ValueError, match="method"): + dq.clean.outlier_mask(x, "bogus") + with pytest.raises(ValueError, match="side"): + dq.clean.outlier_mask(x, "mad", side="above") + # despike: interp bridges the spike, median uses the global median + y = np.array([1.0, 1.0, 50.0, 1.0, 1.0]) + np.testing.assert_allclose(dq.clean.despike(y), np.ones(5)) + np.testing.assert_allclose(dq.clean.despike(y, replace="median"), np.ones(5)) + assert np.isnan(dq.clean.despike(y, replace="nan")[2]) + with pytest.raises(ValueError, match="replace"): + dq.clean.despike(y, replace="zero") + + +def test_closure_renormalize_discrepancy() -> None: + np.testing.assert_allclose(dq.clean.closure_residual(0.3, 0.5, 0.2), 0.0, atol=1e-15) + np.testing.assert_allclose(float(dq.clean.closure_residual(0.6, target=0.5)), 0.1) + with pytest.raises(ValueError, match="at least one"): + dq.clean.closure_residual() + closed = dq.clean.renormalize([[2.0, 2.0], [1.0, 3.0]]) + np.testing.assert_allclose(closed.sum(axis=-1), [1.0, 1.0]) + # deviation-indicator anchor (src2016_02 Eq. 6): 2*0.2/2.2 + np.testing.assert_allclose(dq.clean.relative_discrepancy(1.0, 1.2), 2.0 * 0.2 / 2.2) + + +# --- scale ---------------------------------------------------------------------- + + +def test_scale() -> None: + x = np.linspace(0.0, 10.0, 101) + # explicit endpoints: affine two-point map hits the reference endpoints + out = dq.scale.normalize_to_reference(x, 100.0, 200.0, in_lo=0.0, in_hi=10.0) + np.testing.assert_allclose([out[0], out[-1]], [100.0, 200.0]) + # pct=None uses min/max, so the full range maps exactly + out = dq.scale.normalize_to_reference(x, 0.0, 1.0, pct=None) + np.testing.assert_allclose([out.min(), out.max()], [0.0, 1.0]) + # default 5/95 percentiles: the p5/p95 samples land on the endpoints + out = dq.scale.normalize_to_reference(x, 0.0, 1.0) + np.testing.assert_allclose(out[5], 0.0, atol=1e-12) + np.testing.assert_allclose(out[95], 1.0, atol=1e-12) + y = dq.scale.match_moments(np.random.default_rng(0).normal(5.0, 2.0, 500), 75.0, 20.0) + np.testing.assert_allclose([y.mean(), y.std()], [75.0, 20.0]) + + +# --- filt ----------------------------------------------------------------------- + + +def test_smooth() -> None: + x = np.ones(30) + x[15] = 31.0 + # boxcar: mean-preserving pulse spread; window 1 is a no-op copy + np.testing.assert_allclose(dq.filt.smooth(x, 3)[14:17], [11.0, 11.0, 11.0]) + assert np.array_equal(dq.filt.smooth(x, 1), x) + # median kills the single spike outright + np.testing.assert_allclose(dq.filt.smooth(x, 3, "median"), np.ones(30)) + # gaussian: unit-normalized kernel preserves a constant interior + g = dq.filt.smooth(np.full(50, 7.0), kind="gaussian", sigma=2.0) + np.testing.assert_allclose(g[10:40], 7.0) + with pytest.raises(ValueError, match="sigma"): + dq.filt.smooth(x, kind="gaussian") + with pytest.raises(ValueError, match="odd"): + dq.filt.smooth(x, 4, "median") + with pytest.raises(ValueError, match="kind"): + dq.filt.smooth(x, 5, "hann") + + +def test_moving_stat() -> None: + x = np.arange(10.0) + assert np.array_equal(dq.filt.moving_stat(np.full(20, 3.0), 5, "std"), np.zeros(20)) + # centered window mean of a ramp is the ramp (away from edges) + np.testing.assert_allclose(dq.filt.moving_stat(x, 3, "mean")[1:-1], x[1:-1]) + # causal window lags: trailing mean at i=9 is mean(7,8,9)=8 + np.testing.assert_allclose(dq.filt.moving_stat(x, 3, "mean", center=False)[-1], 8.0) + # cv guards a zero mean with 0 + assert dq.filt.moving_stat(np.zeros(5), 3, "cv")[2] == 0.0 + np.testing.assert_allclose( + dq.filt.moving_stat(x, 3, "var")[5], dq.filt.moving_stat(x, 3, "std")[5] ** 2 + ) + assert dq.filt.moving_stat(np.array([0.0, 4.0, 0.0]), 3, "mad")[1] > 0 + with pytest.raises(ValueError, match="stat"): + dq.filt.moving_stat(x, 3, "median") + + +def test_tool_response_and_median2d() -> None: + # unit-area Gaussian kernel preserves a constant log interior + out = dq.filt.tool_response(np.full(100, 2.5), dz=0.1, fwhm=0.5) + np.testing.assert_allclose(out[20:80], 2.5) + # a step gets smoothed monotonically across the edge + step = np.repeat([0.0, 1.0], 50) + sm = dq.filt.tool_response(step, dz=0.1, fwhm=0.8) + assert 0.0 < sm[49] < sm[50] < sm[51] < 1.0 + # median filter removes salt-and-pepper without moving the background + img = np.full((8, 8), 5.0) + img[3, 4] = 500.0 + np.testing.assert_allclose(dq.filt.median_filter2d(img), np.full((8, 8), 5.0)) + with pytest.raises(ValueError, match="2-D"): + dq.filt.median_filter2d(np.zeros(5)) + with pytest.raises(ValueError, match="odd"): + dq.filt.median_filter2d(img, size=4) + + +def test_window_features_and_boundaries() -> None: + x = np.arange(20.0) + feats = dq.filt.window_features(x, 5) + assert feats.shape == (20, 5) + np.testing.assert_allclose(feats[10], [8.0, 9.0, 10.0, 11.0, 12.0]) + np.testing.assert_allclose(feats[0], [0.0, 0.0, 0.0, 1.0, 2.0]) # edge-padded + both = dq.filt.window_features(np.column_stack([x, 2.0 * x]), 5) + assert both.shape == (20, 10) + with pytest.raises(ValueError, match="odd"): + dq.filt.window_features(x, 4) + # two clean steps -> boundaries at (or next to) the true edges + curve = np.concatenate([np.full(40, 10.0), np.full(30, 80.0), np.full(50, 30.0)]) + curve += np.random.default_rng(1).normal(0.0, 0.5, curve.size) + idx = dq.filt.detect_bed_boundaries(curve, 11) + assert any(abs(i - 40) <= 2 for i in idx) and any(abs(i - 70) <= 2 for i in idx) + + +# --- signal --------------------------------------------------------------------- + + +def test_signal() -> None: + sig = np.full(100, 2.0) + # P_sig=4, sigma=0.2 -> 10*log10(4/0.04) = 20 dB + np.testing.assert_allclose(dq.signal.snr_db(sig, noise_std=0.2), 20.0) + assert dq.signal.snr_db(sig, noise_std=0.0) == float("inf") + np.testing.assert_allclose(dq.signal.snr_db(sig, noise=np.full(100, 0.2)), 20.0) + with pytest.raises(ValueError, match="exactly one"): + dq.signal.snr_db(sig) + with pytest.raises(ValueError, match="exactly one"): + dq.signal.snr_db(sig, noise=sig, noise_std=1.0) + np.testing.assert_allclose(dq.signal.stack_gain(16), 4.0) + # block means over n=3; the trailing partial block (7.0) is dropped + out = dq.signal.block_stack(np.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]), 3) + np.testing.assert_allclose(out, [2.0, 5.0]) + two = dq.signal.block_stack(np.arange(12.0).reshape(2, 6), 3, axis=-1) + np.testing.assert_allclose(two, [[1.0, 4.0], [7.0, 10.0]]) + with pytest.raises(ValueError, match=">= 1"): + dq.signal.block_stack(sig, 0) + # deterministic default rng; sigma scales with the mean absolute level + x = np.full(1000, 10.0) + n1 = dq.signal.add_gaussian_noise(x, 0.05) + n2 = dq.signal.add_gaussian_noise(x, 0.05) + assert np.array_equal(n1, n2) + np.testing.assert_allclose(np.std(n1 - x), 0.5, rtol=0.1) + + +# --- synth ---------------------------------------------------------------------- + + +def test_synth_logs() -> None: + curve, bounds = dq.synth.blocky_log(300, n_beds=6, rng=2) + assert curve.shape == (300,) and bounds.shape == (5,) + assert np.all(np.diff(bounds) > 0) + with pytest.raises(ValueError, match="n_beds"): + dq.synth.blocky_log(300, n_beds=0) + logs, bp, fac = dq.synth.log_suite() + assert set(logs) == set(dq.synth.LOG_SUITE_CURVES) and len(bp) == 3 + assert logs["GR"].shape == (600,) and set(np.unique(fac)) <= {0, 1, 2, 3} + # shale (facies 1) has the highest mean GR of the default properties + assert logs["GR"][fac == 1].mean() > logs["GR"][fac == 0].mean() + with pytest.raises(ValueError, match="n_facies"): + dq.synth.log_suite(n_facies=9) + ref, other = dq.synth.shifted_pair(200, 7) + np.testing.assert_allclose(other[:-7], ref[7:]) + ref, other = dq.synth.shifted_pair(200, -5) + np.testing.assert_allclose(other[5:], ref[:-5]) + + +def test_synth_spectra_and_images() -> None: + t2 = np.logspace(-1.0, 4.0, 200) + spec = dq.synth.gaussian_mixture_spectrum(t2, [1.0, 100.0], [0.6, 0.4], [0.3, 0.3]) + # peaks land at the component centres on the log axis + assert abs(np.log10(t2[np.argmax(spec)]) - 0.0) < 0.05 + np.testing.assert_allclose(spec.max(), 0.6, rtol=1e-3) + noisy = dq.synth.gaussian_mixture_spectrum(t2, [1.0], [0.5], [0.3], noise=0.01, rng=0) + assert noisy.min() >= 0.0 + with pytest.raises(ValueError): + dq.synth.gaussian_mixture_spectrum(t2, [1.0, 2.0], [0.5], [0.3]) + cube = dq.synth.sphere_pack_volume((24, 24, 24), (3, 5), 0.4, rng=0) + phi = 1.0 - cube.mean() + assert cube.dtype == np.uint8 and phi <= 0.4 + img, mask = dq.synth.disk_image(64, 20, rng=4) + assert img.shape == (64, 64) and mask.dtype == bool + assert 0.0 <= img.min() and img.max() <= 1.0 + # disks are darker than the background + assert img[mask].mean() < img[~mask].mean() + + +# --- io ------------------------------------------------------------------------- + + +def test_wellbore_data(tmp_path: Path) -> None: + w = dq.WellboreData("Well-A", uwi="00/01") + w.add_channel("GR", np.linspace(10.0, 90.0, 8), "gAPI", ["MD"]) + w.add_channel("UDAR", np.ones((4, 3, 2)), "ohm.m", ["MD", "azimuth", "DOI"]) + with pytest.raises(ValueError, match="axes count"): + w.add_channel("BAD", np.ones((2, 2)), "-", ["MD"]) + s = w.to_json(str(tmp_path / "well.json")) + assert (tmp_path / "well.json").read_text() == s + w2 = dq.WellboreData.from_json(s) + assert w2.meta == w.meta + assert w2.channels["UDAR"]["shape"] == (4, 3, 2) + np.testing.assert_array_equal(w2.channels["GR"]["data"], w.channels["GR"]["data"])