diff --git a/src2023_04/article10_induction_deconvolution.py b/src2023_04/article10_induction_deconvolution.py index 34639d0..9f0b7fd 100644 --- a/src2023_04/article10_induction_deconvolution.py +++ b/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) diff --git a/src2023_04/article11_induction_convolution.py b/src2023_04/article11_induction_convolution.py index a433e16..3405a98 100644 --- a/src2023_04/article11_induction_convolution.py +++ b/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) diff --git a/src2023_06/article9_xray_invasion.py b/src2023_06/article9_xray_invasion.py index e69aefe..c567c9b 100644 --- a/src2023_06/article9_xray_invasion.py +++ b/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) # --------------------------------------------------------------------------- diff --git a/src2023_08/article5_iterative_resistivity.py b/src2023_08/article5_iterative_resistivity.py index cefd5ac..066dadd 100644 --- a/src2023_08/article5_iterative_resistivity.py +++ b/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) # --------------------------------------------------------------------------- diff --git a/src2023_10/article_03_okwoli_probe_screening.py b/src2023_10/article_03_okwoli_probe_screening.py index 9208a13..cac2288 100644 --- a/src2023_10/article_03_okwoli_probe_screening.py +++ b/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, diff --git a/src2023_10/article_09_bennis_corelogs_simulation.py b/src2023_10/article_09_bennis_corelogs_simulation.py index 0f26f7e..6e192e2 100644 --- a/src2023_10/article_09_bennis_corelogs_simulation.py +++ b/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) # --------------------------------------------------------------------------- diff --git a/src2023_12/bradley_wellbore_format.py b/src2023_12/bradley_wellbore_format.py index 41c3414..4349d64 100644 --- a/src2023_12/bradley_wellbore_format.py +++ b/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(): diff --git a/src2024_04/chen_sem_pore_segmentation.py b/src2024_04/chen_sem_pore_segmentation.py index 3ea87e4..e6890ef 100644 --- a/src2024_04/chen_sem_pore_segmentation.py +++ b/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): diff --git a/src2024_08/magnetic_permeability.py b/src2024_08/magnetic_permeability.py index 1f294f7..674577a 100644 --- a/src2024_08/magnetic_permeability.py +++ b/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, diff --git a/src2025_08/anomaly_detection_vmd.py b/src2025_08/anomaly_detection_vmd.py index 56b28a4..e488cf7 100644 --- a/src2025_08/anomaly_detection_vmd.py +++ b/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) # --------------------------------------------------------------------------- diff --git a/src2025_08/fiber_optics_sensing.py b/src2025_08/fiber_optics_sensing.py index 909d0e1..537346d 100644 --- a/src2025_08/fiber_optics_sensing.py +++ b/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) # --------------------------------------------------------------------------- diff --git a/src2026_04/a04_wireline_anomaly_diagnosis.py b/src2026_04/a04_wireline_anomaly_diagnosis.py index 2e60f35..6f88ded 100644 --- a/src2026_04/a04_wireline_anomaly_diagnosis.py +++ b/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: diff --git a/src2026_06/a03_pore_size_fluid_movement.py b/src2026_06/a03_pore_size_fluid_movement.py index 60773c0..f9397af 100644 --- a/src2026_06/a03_pore_size_fluid_movement.py +++ b/src2026_06/a03_pore_size_fluid_movement.py @@ -105,11 +105,13 @@ class GaussianMode: def gaussian_t2_model(log_t2: np.ndarray, modes: Sequence[GaussianMode]) -> np.ndarray: """Sum of Gaussian pore-body modes evaluated on a log10(T2) axis.""" - log_t2 = np.asarray(log_t2, dtype=float) - out = np.zeros_like(log_t2) - for m in modes: - out += m.amplitude * np.exp(-0.5 * ((log_t2 - m.center) / m.width) ** 2) - return out + return petrolib.data_qc_io.synth.gaussian_mixture_spectrum( + log_t2, + [m.center for m in modes], + [m.amplitude for m in modes], + [m.width for m in modes], + log_axis=False, + ) def partition_porosity(t2_ms: np.ndarray, amplitude: np.ndarray,