diff --git a/src2020_10/article3_cavein_clustering_detection.py b/src2020_10/article3_cavein_clustering_detection.py
index f6066c2..09c36f4 100644
--- a/src2020_10/article3_cavein_clustering_detection.py
+++ b/src2020_10/article3_cavein_clustering_detection.py
@@ -39,15 +39,7 @@
 
 def rolling_cv(rhob, wr):
     """Rolling coefficient of variation of bulk density over +/- wr  (Eq. 5)."""
-    rhob = np.asarray(rhob, float)
-    n = len(rhob)
-    cv = np.zeros(n)
-    for i in range(n):
-        lo, hi = max(0, i - wr), min(n, i + wr + 1)
-        w = rhob[lo:hi]
-        m = w.mean()
-        cv[i] = w.std() / m if m != 0 else 0.0
-    return cv
+    return petrolib.data_qc_io.filt.moving_stat(rhob, 2 * wr + 1, "cv")
 
 
 def caliper_excess(caliper, bit_size):
diff --git a/src2021_04/article5_nmr_fluid_contacts.py b/src2021_04/article5_nmr_fluid_contacts.py
index 309966b..245ce0e 100644
--- a/src2021_04/article5_nmr_fluid_contacts.py
+++ b/src2021_04/article5_nmr_fluid_contacts.py
@@ -88,7 +88,7 @@ def classify_fluid(D, gas_D=1e-8, oil_D=5e-10):
 
 def stacking_snr(base_snr, n_stacks):
     """Station stacking improves SNR as sqrt(number of stacked acquisitions)."""
-    return base_snr * np.sqrt(n_stacks)
+    return float(base_snr * petrolib.data_qc_io.signal.stack_gain(n_stacks))
 
 
 # ---------------------------------------------- tests --------------
diff --git a/src2021_12/article01_data_quality_ml.py b/src2021_12/article01_data_quality_ml.py
index d2d1187..aeae8b7 100644
--- a/src2021_12/article01_data_quality_ml.py
+++ b/src2021_12/article01_data_quality_ml.py
@@ -51,7 +51,7 @@ def zscore(x):
 
 def zscore_outliers(x, threshold=3.0):
     """Boolean outlier mask using |z| > threshold (default 3 -> 99.7%)."""
-    return np.abs(zscore(x)) > threshold
+    return petrolib.data_qc_io.clean.outlier_mask(x, "zscore", threshold=threshold)
 
 
 # ---------------------------------------------- Eq. 2: IQR --------------
@@ -64,9 +64,7 @@ def iqr_bounds(x, k=1.5):
 
 
 def iqr_outliers(x, k=1.5):
-    lo, hi = iqr_bounds(x, k)
-    x = np.asarray(x, float)
-    return (x < lo) | (x > hi)
+    return petrolib.data_qc_io.clean.outlier_mask(x, "iqr", k=k)
 
 
 # ---------------------------------------------- Eqs. 3-4: normalize -----
@@ -82,8 +80,8 @@ def normalize_reference(v, R_min, R_max, W_min, W_max):
     V_norm = R_min + (R_max - R_min) * (V - W_min)/(W_max - W_min).
     R_* are reference-well 5th/95th percentiles, W_* the target well's.
     """
-    v = np.asarray(v, float)
-    return R_min + (R_max - R_min) * (v - W_min) / (W_max - W_min)
+    return petrolib.data_qc_io.scale.normalize_to_reference(
+        v, R_min, R_max, in_lo=W_min, in_hi=W_max)
 
 
 # ---------------------------------------------- Eqs. 5-6: P / R ---------
@@ -114,10 +112,7 @@ def rmse(y_true, y_pred):
 
 def add_gaussian_noise(x, sigma_fraction, rng=None):
     """Add N(0, sigma^2) noise with sigma = fraction * mean(|x|) (Eqs. 9-10)."""
-    x = np.asarray(x, float)
-    rng = rng or np.random.default_rng(0)
-    sigma = sigma_fraction * np.abs(x).mean()
-    return x + rng.normal(0.0, sigma, size=x.shape)
+    return petrolib.data_qc_io.signal.add_gaussian_noise(x, sigma_fraction, rng)
 
 
 # ---------------------------------------------- Eq. 11: Pearson ---------
@@ -132,10 +127,7 @@ def pearson(x, y):
 
 def sentinels_to_nan(x, sentinels=SENTINELS):
     """Replace LAS/DLIS no-data sentinels with NaN."""
-    x = np.asarray(x, float).copy()
-    for s in sentinels:
-        x[np.isclose(x, s)] = np.nan
-    return x
+    return petrolib.data_qc_io.clean.sentinels_to_nan(x, sentinels)
 
 
 # ---------------------------------------------- tests --------------
diff --git a/src2022_08/article5_cpa_fcm_logfacies.py b/src2022_08/article5_cpa_fcm_logfacies.py
index 3561681..8f031ef 100644
--- a/src2022_08/article5_cpa_fcm_logfacies.py
+++ b/src2022_08/article5_cpa_fcm_logfacies.py
@@ -25,6 +25,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
+
 
 # -------------------------------------------- CPA (Eqs. 1-5) --------------
 
@@ -135,26 +142,10 @@ def fcm(X, c=5, m=2.0, n_iter=200, tol=1e-5, seed=0):
 
 def synth_log_suite(n=600, n_facies=4, seed=0):
     """Five-log suite (GR, RD, DEN, AT, NP) with `n_facies` blocks."""
-    rng = np.random.default_rng(seed)
-    bp_true = sorted(rng.choice(np.arange(40, n - 40), n_facies - 1,
-                                replace=False))
-    bp_full = [0] + list(bp_true) + [n]
-    facies_props = np.array([
-        [ 22.0,  80.0, 2.45, 75.0, 0.08],   # facies 0: granitic, low GR
-        [ 75.0,  10.0, 2.20, 95.0, 0.30],   # facies 1: shale
-        [ 35.0, 200.0, 2.60, 65.0, 0.05],   # facies 2: fractured granite
-        [ 50.0,  30.0, 2.40, 80.0, 0.15],   # facies 3: weathered
-    ])
-    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:])):
-        fac = k % n_facies
-        facies_labels[a:b] = fac
-        for j in range(5):
-            mu = facies_props[fac, j]
-            sigma = 0.05 * abs(mu) + 1e-3
-            out[a:b, j] = rng.normal(mu, sigma, b - a)
-    return out, np.array(bp_true), facies_labels
+    logs, bp_true, facies_labels = petrolib.data_qc_io.synth.log_suite(
+        n, n_facies, rng=np.random.default_rng(seed))
+    out = np.column_stack([logs[c] for c in petrolib.data_qc_io.synth.LOG_SUITE_CURVES])
+    return out, bp_true, facies_labels
 
 
 # -------------------------------------------- tests ----------------------