Adopt petrolib.data_qc_io in 2023-2026 articles (PR-B3, final chunk of the merge plan)

src2023_04/article10_induction_deconvolution.py

@@ -14,6 +14,13 @@
 """
 
 import numpy as np
+
+try:
+    import petrolib
+except ImportError:  # bare clone, not installed
+    import sys, pathlib
+    sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
+    import petrolib
 import xgboost as xgb
 from scipy.ndimage import gaussian_filter1d
 
@@ -59,10 +66,7 @@ def linear_deconvolution(R_app, weights):
 
 def fit_linear_weights(R_app, R_model, window=21):
     """Fit linear regression weights minimising log-RMSE."""
-    pad = window // 2
-    log_app = np.log(R_app)
-    log_app_pad = np.pad(log_app, pad, mode="edge")
-    rows = np.lib.stride_tricks.sliding_window_view(log_app_pad, window)
+    rows = petrolib.data_qc_io.filt.window_features(np.log(R_app), window)
     coef, *_ = np.linalg.lstsq(rows, np.log(R_model), rcond=None)
     return coef
 
@@ -79,10 +83,8 @@ def __init__(self, window=21, n_estimators=300, max_depth=5, seed=0):
                                       verbosity=0)
 
     def _features(self, R_app):
-        pad = self.window // 2
         log_app = np.log(np.maximum(R_app, 1e-9))
-        log_pad = np.pad(log_app, pad, mode="edge")
-        return np.lib.stride_tricks.sliding_window_view(log_pad, self.window).copy()
+        return petrolib.data_qc_io.filt.window_features(log_app, self.window)
 
     def fit(self, R_app, R_model):
         X = self._features(R_app)

src2023_04/article11_induction_convolution.py

@@ -12,6 +12,13 @@
 """
 
 import numpy as np
+
+try:
+    import petrolib
+except ImportError:  # bare clone, not installed
+    import sys, pathlib
+    sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
+    import petrolib
 import xgboost as xgb
 
 from article10_induction_deconvolution import (
@@ -24,10 +31,8 @@
 # ----------------------- linear convolution baseline (Eq. 2) ---
 
 def fit_linear_conv_weights(R_model, R_app, window=101):
-    pad = window // 2
     log_model = np.log(np.maximum(R_model, 1e-9))
-    log_pad = np.pad(log_model, pad, mode="edge")
-    rows = np.lib.stride_tricks.sliding_window_view(log_pad, window)
+    rows = petrolib.data_qc_io.filt.window_features(log_model, window)
     coef, *_ = np.linalg.lstsq(rows, np.log(np.maximum(R_app, 1e-9)), rcond=None)
     return coef
 
@@ -52,10 +57,8 @@ def __init__(self, window=101, n_estimators=300, max_depth=5, seed=0):
                                       verbosity=0)
 
     def _features(self, R_model):
-        pad = self.window // 2
         log_m = np.log(np.maximum(R_model, 1e-9))
-        log_pad = np.pad(log_m, pad, mode="edge")
-        return np.lib.stride_tricks.sliding_window_view(log_pad, self.window).copy()
+        return petrolib.data_qc_io.filt.window_features(log_m, self.window)
 
     def fit(self, R_model, R_app):
         X = self._features(R_model)

src2023_06/article9_xray_invasion.py

@@ -36,6 +36,13 @@
 
 import numpy as np
 
+try:
+    import petrolib
+except ImportError:  # bare clone, not installed
+    import sys, pathlib
+    sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
+    import petrolib
+
 
 # ---------------------------------------------------------------------------
 # Beer-Lambert
@@ -57,15 +64,7 @@ def attenuation_map(I_now: np.ndarray, I_dry: np.ndarray) -> np.ndarray:
 # ---------------------------------------------------------------------------
 def median_filter_3x3(image: np.ndarray) -> np.ndarray:
     """Pure-numpy 3x3 median filter (with edge replication)."""
-    pad = np.pad(image, 1, mode="edge")
-    H, W = image.shape
-    stack = np.empty((9, H, W))
-    k = 0
-    for dy in (0, 1, 2):
-        for dx in (0, 1, 2):
-            stack[k] = pad[dy:dy + H, dx:dx + W]
-            k += 1
-    return np.median(stack, axis=0)
+    return petrolib.data_qc_io.filt.median_filter2d(image, size=3)
 
 
 # ---------------------------------------------------------------------------

src2023_08/article5_iterative_resistivity.py

@@ -39,6 +39,13 @@
 
 import numpy as np
 
+try:
+    import petrolib
+except ImportError:  # bare clone, not installed
+    import sys, pathlib
+    sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
+    import petrolib
+
 
 # ---------------------------------------------------------------------------
 # 1. Bed-boundary detection (variance + first-derivative sliding window)
@@ -50,23 +57,8 @@ def detect_bed_boundaries(log_curve, depth, win=11, threshold=None):
     sliding window computes (a) the first derivative magnitude and
     (b) the local variance of the log; their product is thresholded.
     """
-    log_curve = np.asarray(log_curve, dtype=float)
-    n = len(log_curve)
-    half = win // 2
-    deriv = np.gradient(log_curve)
-    variance = np.zeros(n)
-    for i in range(half, n - half):
-        variance[i] = np.var(log_curve[i - half:i + half + 1])
-    metric = np.abs(deriv) * np.sqrt(variance)
-    if threshold is None:
-        threshold = np.mean(metric) + 2.0 * np.std(metric)
-    # Local maxima of metric above threshold
-    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)
+    return petrolib.data_qc_io.filt.detect_bed_boundaries(
+        log_curve, window=win, threshold=threshold)
 
 
 # ---------------------------------------------------------------------------

src2023_10/article_03_okwoli_probe_screening.py

@@ -42,6 +42,13 @@
 
 import numpy as np
 
+try:
+    import petrolib
+except ImportError:  # bare clone, not installed
+    import sys, pathlib
+    sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
+    import petrolib
+
 
 # ---------------------------------------------------------------------------
 # Synthetic core
@@ -99,10 +106,7 @@ def make_synthetic_core(length_m: float = 5.0,
 # ---------------------------------------------------------------------------
 def boxcar(x: np.ndarray, window: int) -> np.ndarray:
     """Centered moving average."""
-    if window <= 1:
-        return x.copy()
-    k = np.ones(window) / window
-    return np.convolve(x, k, mode="same")
+    return petrolib.data_qc_io.filt.smooth(x, window, "boxcar")
 
 
 def upscale_to(probe: CoreProbes, sample_dz_mm: float,

src2023_10/article_09_bennis_corelogs_simulation.py

@@ -115,12 +115,7 @@ def synthetic_logs(s: CoreSample) -> dict:
 # ---------------------------------------------------------------------------
 def vertical_response(x: np.ndarray, dz_m: float, fwhm_m: float) -> np.ndarray:
     """Convolve x with a Gaussian whose FWHM matches the tool aperture."""
-    sigma = fwhm_m / (2.0 * np.sqrt(2.0 * np.log(2.0)))
-    n = max(3, int(round(6 * sigma / dz_m)) | 1)
-    t = (np.arange(n) - n // 2) * dz_m
-    g = np.exp(-0.5 * (t / sigma) ** 2)
-    g /= g.sum()
-    return np.convolve(x, g, mode="same")
+    return petrolib.data_qc_io.filt.tool_response(x, dz_m, fwhm_m)
 
 
 # ---------------------------------------------------------------------------

src2023_12/bradley_wellbore_format.py

@@ -10,37 +10,16 @@
 import json
 import numpy as np
 
+try:
+    import petrolib
+except ImportError:  # bare clone, not installed
+    import sys, pathlib
+    sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
+    import petrolib
 
-class WellboreData:
-    def __init__(self, well_name, uwi=None):
-        self.meta = {"well_name": well_name, "uwi": uwi, "version": "1.0"}
-        self.channels = {}  # name -> dict(data, units, axes, dims)
-
-    def add_channel(self, name, data, units, axes):
-        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):
-        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=None):
-        s = json.dumps(self.to_dict(), indent=2)
-        if path: open(path, "w").write(s)
-        return s
-
-    @classmethod
-    def from_dict(cls, d):
-        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
+
+class WellboreData(petrolib.data_qc_io.WellboreData):
+    """JSON wellbore container -- canonical implementation in petrolib."""
 
 
 def test_all():

src2024_04/chen_sem_pore_segmentation.py

@@ -15,6 +15,13 @@
 """
 import numpy as np
 
+try:
+    import petrolib
+except ImportError:  # bare clone, not installed
+    import sys, pathlib
+    sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
+    import petrolib
+
 try:
     from scipy import ndimage as ndi
 except Exception as e:
@@ -23,16 +30,7 @@
 
 def synthetic_sem_image(size=128, n_pores=40, rng=None):
     rng = rng or np.random.default_rng(4)
-    img = 0.6 + 0.05 * rng.standard_normal((size, size))
-    mask = np.zeros((size, size), dtype=bool)
-    yy, xx = np.mgrid[:size, :size]
-    for _ in range(n_pores):
-        cy, cx = rng.integers(0, size, size=2)
-        r = rng.integers(2, 8)
-        m = (yy - cy) ** 2 + (xx - cx) ** 2 <= r * r
-        mask |= m
-    img[mask] = 0.15 + 0.05 * rng.standard_normal(mask.sum())
-    return np.clip(img, 0, 1), mask
+    return petrolib.data_qc_io.synth.disk_image(size, n_pores, rng=rng)
 
 
 def porenet_segment(img, threshold=0.35):

src2024_08/magnetic_permeability.py

@@ -15,6 +15,13 @@
 """
 
 import numpy as np
+
+try:
+    import petrolib
+except ImportError:  # bare clone, not installed
+    import sys, pathlib
+    sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
+    import petrolib
 from dataclasses import dataclass
 from typing import Tuple, Optional
 from scipy.ndimage import gaussian_filter1d
@@ -62,14 +69,8 @@ def remove_ferromagnetic_contaminants(susceptibility: np.ndarray,
     very high susceptibilities that saturate in high fields, but
     residual signals may remain.
     """
-    median_val = np.median(susceptibility)
-    mad = np.median(np.abs(susceptibility - median_val))
-    upper_bound = median_val + threshold_factor * mad * 1.4826
-    cleaned = susceptibility.copy()
-    outliers = cleaned > upper_bound
-    if outliers.any():
-        cleaned[outliers] = median_val
-    return cleaned
+    return petrolib.data_qc_io.clean.despike(
+        susceptibility, threshold=threshold_factor, side="upper", replace="median")
 
 
 def estimate_clay_volume_illite(susceptibility_total: np.ndarray,

src2025_08/anomaly_detection_vmd.py

@@ -18,6 +18,13 @@
 """
 
 import numpy as np
+
+try:
+    import petrolib
+except ImportError:  # bare clone, not installed
+    import sys, pathlib
+    sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
+    import petrolib
 from typing import List, Tuple, Optional
 from dataclasses import dataclass
 
@@ -262,10 +269,8 @@ def compute_snr_improvement(
     original: np.ndarray, mode_signal: np.ndarray,
 ) -> float:
     """Estimate SNR improvement (dB) of extracted mode vs. original."""
-    noise_est = original - mode_signal
-    sig_power = np.mean(mode_signal ** 2)
-    noise_power = np.mean(noise_est ** 2) + 1e-30
-    return 10.0 * np.log10(sig_power / noise_power)
+    return petrolib.data_qc_io.signal.snr_db(
+        mode_signal, noise=original - mode_signal, guard=1e-30)
 
 
 # ---------------------------------------------------------------------------

src2025_08/fiber_optics_sensing.py

@@ -16,6 +16,13 @@
 """
 
 import numpy as np
+
+try:
+    import petrolib
+except ImportError:  # bare clone, not installed
+    import sys, pathlib
+    sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
+    import petrolib
 from typing import List, Tuple, Optional
 from dataclasses import dataclass
 
@@ -124,21 +131,12 @@ def temporal_stack(das: np.ndarray, n_stack: int = 10) -> np.ndarray:
 
     The SNR improves approximately as sqrt(n_stack).
     """
-    n_depth, n_time = das.shape
-    n_out = n_time // n_stack
-    stacked = np.zeros((n_depth, n_out))
-    for i in range(n_out):
-        stacked[:, i] = das[:, i * n_stack:(i + 1) * n_stack].mean(axis=1)
-    return stacked
+    return petrolib.data_qc_io.signal.block_stack(das, n_stack, axis=-1)
 
 
 def compute_snr(signal: np.ndarray, noise_std: float) -> float:
     """Compute signal-to-noise ratio in dB."""
-    sig_power = np.mean(signal ** 2)
-    noise_power = noise_std ** 2
-    if noise_power == 0:
-        return float('inf')
-    return 10.0 * np.log10(sig_power / noise_power)
+    return petrolib.data_qc_io.signal.snr_db(signal, noise_std=noise_std)
 
 
 # ---------------------------------------------------------------------------

src2026_04/a04_wireline_anomaly_diagnosis.py

@@ -119,29 +119,17 @@ def winch_acceleration(tension_variation: np.ndarray,
 
 def moving_average(x: np.ndarray, window: int) -> np.ndarray:
     """Moving average MA_t over past N points (Eq. 13)."""
-    kernel = np.ones(window) / window
-    return np.convolve(x, kernel, mode='same')
+    return petrolib.data_qc_io.filt.smooth(x, window, "boxcar")
 
 
 def moving_msd(x: np.ndarray, window: int) -> np.ndarray:
     """Moving mean-squared deviation MSD_t (Eq. 14)."""
-    n = len(x)
-    out = np.zeros(n)
-    for i in range(n):
-        lo = max(0, i - window + 1)
-        seg = x[lo:i+1]
-        out[i] = np.mean((seg - seg.mean())**2)
-    return out
+    return petrolib.data_qc_io.filt.moving_stat(x, window, "var", center=False)
 
 
 def moving_std(x: np.ndarray, window: int) -> np.ndarray:
     """Moving standard deviation σ_t (Eq. 15)."""
-    n = len(x)
-    out = np.zeros(n)
-    for i in range(n):
-        lo = max(0, i - window + 1)
-        out[i] = np.std(x[lo:i+1])
-    return out
+    return petrolib.data_qc_io.filt.moving_stat(x, window, "std", center=False)
 
 
 def moving_abs_deviation(x: np.ndarray, ma: np.ndarray) -> np.ndarray: