Examples

The examples/ directory in the source tree contains four runnable scripts. Each is self-contained and uses the synthetic-data generator so it can run on any machine without external files.

quick_start.py

A three-line tour of the package: load synthetic data, auto-detect events, run Bourdet / Horner / reservoir-parameter analysis on the longest buildup.

"""
quick_start.py — three-line tour of welltest-pta.

This example uses synthetic data so it can run without any external file.
Replace the synthetic load with WellTest.from_file("yourfile.txt") for real
gauge data.
"""

from welltest_pta import WellTest
from welltest_pta.utils.synthetic import generate_synthetic_dst


# ── 1. Load (here: synthetic; in practice: from a gauge file) ────────────
df = generate_synthetic_dst(n_samples=10_000, sample_period_s=4.0)
wt = WellTest.from_dataframe(df)

# ── 2. Inspect the event catalogue ───────────────────────────────────────
wt.print_summary()

# ── 3. Per-event analysis on the longest buildup ─────────────────────────
bu = wt.events.longest_buildup
bu.print()

print("\n>>> Bourdet log-derivative")
dt, deriv = bu.bourdet(L=0.2)
print(f"    {len(dt)} derivative points covering Δt = "
      f"[{dt[0]:.4f}, {dt[-1]:.2f}] hr")

print("\n>>> Horner extrapolation")
h = bu.horner()
print(f"    P*    = {h['p_star']:.1f} psi")
print(f"    slope = {h['slope_m']:.2f} psi/cycle")
print(f"    R²    = {h['r2']:.4f}")

print("\n>>> Reservoir parameters (oilfield units)")
params = bu.reservoir_params(
    q=850, mu=0.45, B=1.18,
    h=18, phi=0.12, ct=1.2e-5, rw=0.108,
    method="horner",
)
print(f"    k    = {params['k']:.3f} mD")
print(f"    kh   = {params['kh']:.1f} mD·ft")
print(f"    skin = {params['skin']:+.3f}")

print("\n>>> Flow-regime identification")
for reg in bu.flow_regimes():
    print(f"    {reg['regime']:<22s} slope={reg['slope_mean']:+.2f}  "
          f"Δt = [{reg['dt_start']:.3f}, {reg['dt_end']:.3f}] hr")

manual_split.py

Demonstrates the workflow when the auto-detector’s CV score is marginal. After auto-detection, the user passes a list of (type, t_start, t_end) tuples to split_manual(), which rebuilds the EventCollection from explicit timestamps.

"""
manual_split.py — manual override of the auto-detector.

Demonstrates the workflow when the cross-validation score is marginal:
1. Run auto-detect first (required — sets up p_smooth / elapsed_hr columns).
2. Inspect the catalogue and decide which events you want to override.
3. Pass a list of (type, t_start, t_end) tuples to wt.split_manual().
"""

from welltest_pta import WellTest, EventDetectorConfig
from welltest_pta.utils.synthetic import generate_synthetic_dst


# ── 1. Synthetic test with 4 events (DD-BU-DD-BU) ────────────────────────
df = generate_synthetic_dst(
    n_samples=20_000,
    sample_period_s=3.0,
    sequence=[
        ("DD", 0.5, 3300.0),   # short flow #1
        ("BU", 1.0, 4490.0),   # short BU
        ("DD", 1.0, 3000.0),   # main flow
        ("BU", 8.0, 4495.0),   # extended BU
    ],
)

# ── 2. Auto-detection first ──────────────────────────────────────────────
wt = WellTest.from_dataframe(df, auto_detect=True)
print("[AUTO]")
wt.events.print()

# ── 3. Suppose the auto-detector misclassified BU-1 as part of BU-2.
#     We override with explicit timestamps.
ts = df["timestamp"]

# Find indices roughly bracketing the four events (the synthetic generator
# places them in order with small inter-segment gaps).
idx_dd1_a, idx_dd1_b = 4_500, 5_500
idx_bu1_a, idx_bu1_b = 5_700, 7_500
idx_dd2_a, idx_dd2_b = 7_700, 9_000
idx_bu2_a, idx_bu2_b = 9_100, 18_500

wt.split_manual([
    ("DD", ts.iloc[idx_dd1_a], ts.iloc[idx_dd1_b]),
    ("BU", ts.iloc[idx_bu1_a], ts.iloc[idx_bu1_b]),
    ("DD", ts.iloc[idx_dd2_a], ts.iloc[idx_dd2_b]),
    ("BU", ts.iloc[idx_bu2_a], ts.iloc[idx_bu2_b]),
])

print("\n[MANUAL OVERRIDE]")
wt.events.print()

# ── 4. Per-event analysis is now identical to the auto-detect case ─────
bu_long = wt.events.longest_buildup
print(f"\nLongest BU = {bu_long.event_id}, duration = {bu_long.duration_hr:.2f} hr")
print(f"Preceding tp = {bu_long.preceding_dd_dur_hr:.3f} hr")

h = bu_long.horner()
print(f"Horner: P* = {h['p_star']:.1f} psi, m = {h['slope_m']:.2f}, R² = {h['r2']:.4f}")

deconvolution_demo.py

Multi-event deconvolution on a 6-event synthetic test. The recovered unit-rate response merges all DDs and BUs into one diagnostic plot.

"""
deconvolution_demo.py — multi-event deconvolution merging multiple buildups.

This example shows how deconvolve() takes all DD/BU events from a single
test and merges them into one equivalent unit-rate response — the
"diagnostic master plot" of the entire test.
"""

import matplotlib.pyplot as plt

from welltest_pta import WellTest, deconvolve
from welltest_pta.utils.synthetic import generate_synthetic_dst


# ── 1. Multi-rate synthetic test with three buildups ─────────────────────
df = generate_synthetic_dst(
    n_samples=24_000,
    sample_period_s=4.0,
    sequence=[
        ("DD", 0.4, 3500.0),
        ("BU", 0.8, 4480.0),
        ("DD", 0.8, 3100.0),
        ("BU", 1.5, 4490.0),
        ("DD", 1.0, 2900.0),
        ("BU", 8.0, 4495.0),
    ],
)
wt = WellTest.from_dataframe(df)
wt.print_summary()

# ── 2. Deconvolve ────────────────────────────────────────────────────────
print("\n→ Running vSH04 deconvolution …")
res = deconvolve(
    wt.events,
    default_q=850,                # STB/D for drawdowns
    nu=1e-2,                      # regularisation
    n_response_nodes=60,
    fit_p_initial=True,
)
print(f"   converged = {res.converged}, iters = {res.iterations}")
print(f"   p_initial = {res.p_initial:.1f} psi   ||r|| = {res.residual_norm:.2f} psi")

# ── 3. Plot ──────────────────────────────────────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Left: observed pressure + reconstruction
axes[0].plot(res.obs_time, res.obs_pressure, "o", ms=1.5,
             color="#888888", label="Observed")
axes[0].plot(res.obs_time, res.fit_pressure, "-", lw=1.0,
             color="#1f77b4", label="Reconstructed")
axes[0].set_xlabel("Elapsed time (hr)")
axes[0].set_ylabel("Pressure (psi)")
axes[0].set_title("Pressure: observed vs deconvolution reconstruction")
axes[0].legend()
axes[0].grid(alpha=0.3)

# Right: log-log diagnostic of the recovered unit-rate response
res.plot(ax=axes[1])

plt.tight_layout()
plt.savefig("deconvolution_demo.png", dpi=200, bbox_inches="tight")
print("\n→ Saved deconvolution_demo.png")

full_workflow.py

End-to-end pipeline: load → cross-validate → per-event analytics → reservoir parameters → deconvolution → composite report → bulk export. The recommended starting point for a real analysis.

"""
full_workflow.py — end-to-end pipeline demo.

Shows everything in sequence:
  load → cross-validate → manual override → per-event analytics →
  reservoir parameters → deconvolution → composite report → bulk export.
"""

from pathlib import Path
import warnings
warnings.filterwarnings("ignore")

import matplotlib.pyplot as plt

from welltest_pta import WellTest, deconvolve, EventDetectorConfig
from welltest_pta.utils.synthetic import generate_synthetic_dst


OUTPUT_DIR = Path("welltest_output")
OUTPUT_DIR.mkdir(exist_ok=True)


# ─── 1. Load gauge data ──────────────────────────────────────────────────
print("\n[1/7]  Loading gauge data…")
df = generate_synthetic_dst(
    n_samples=20_000,
    sample_period_s=4.0,
    sequence=[
        ("DD", 0.5, 3300.0),
        ("BU", 1.2, 4485.0),
        ("DD", 1.0, 3000.0),
        ("BU", 8.0, 4495.0),
    ],
)
df.drop(columns="true_event").to_csv(OUTPUT_DIR / "synthetic_input.csv", index=False)


# ─── 2. WellTest construction with cross-validation ─────────────────────
print("\n[2/7]  Auto-detect + cross-validate…")
cfg = EventDetectorConfig(
    hampel_sigma=3.0,
    spike_percentile=95.0,
    min_pta_dp_psi=15.0,
    tail_trim_enabled=True,
)
wt = WellTest.from_dataframe(df, cfg=cfg)
wt.cross_validate(n_bootstrap=4)


# ─── 3. Print catalogue ──────────────────────────────────────────────────
print("\n[3/7]  Event catalogue")
wt.print_summary()


# ─── 4. Per-event analysis ───────────────────────────────────────────────
print("\n[4/7]  Per-event analytics")
for ev in wt.events:
    print(f"\n   {ev}")
    if ev.event_type != "buildup":
        continue
    h = ev.horner()
    print(f"     Horner P*  = {h['p_star']:.1f} psi   m = {h['slope_m']:.2f}  R² = {h['r2']:.4f}")
    m = ev.mdh()
    print(f"     MDH p_1hr  = {m['intercept_p1hr']:.1f} psi   m = {m['slope_m']:.2f}  R² = {m['r2']:.4f}")
    fr = ev.flow_regimes()
    if fr:
        names = ", ".join(r["regime"] for r in fr)
        print(f"     Regimes    : {names}")


# ─── 5. Reservoir parameters on the longest BU ──────────────────────────
print("\n[5/7]  Reservoir parameters (longest BU)")
bu = wt.events.longest_buildup
params = bu.reservoir_params(
    q=850, mu=0.45, B=1.18,
    h=18, phi=0.12, ct=1.2e-5, rw=0.108,
    method="horner",
)
for k, v in params.items():
    if v is None or (isinstance(v, float) and v != v):  # NaN check
        print(f"     {k:<8s} = (n/a)")
    elif isinstance(v, float):
        print(f"     {k:<8s} = {v:+.4f}")
    else:
        print(f"     {k:<8s} = {v}")


# ─── 6. Deconvolution ────────────────────────────────────────────────────
print("\n[6/7]  Multi-event deconvolution")
res = deconvolve(wt.events, default_q=850, nu=1e-2, n_response_nodes=50)
print(f"     converged = {res.converged}, ||r|| = {res.residual_norm:.2f} psi")
res.export(OUTPUT_DIR / "deconvolution.csv")


# ─── 7. Plots + bulk export ──────────────────────────────────────────────
print("\n[7/7]  Plots and bulk export")
fig = wt.plot_composite(out_path=OUTPUT_DIR / "composite_report.pdf")
plt.close(fig)
fig = bu.plot_loglog(); plt.savefig(OUTPUT_DIR / "longest_BU_loglog.png", dpi=180); plt.close(fig)
fig = bu.plot_horner(); plt.savefig(OUTPUT_DIR / "longest_BU_horner.png", dpi=180); plt.close(fig)
fig = res.plot(); plt.savefig(OUTPUT_DIR / "deconvolution.png", dpi=180); plt.close(fig)

paths = wt.export_all(OUTPUT_DIR, prefix="full_workflow", per_event=True)
print(f"\n   Wrote {len(paths)} files to {OUTPUT_DIR}/")
print("\n✓ Pipeline complete.")

Running the examples

After installing the package in development mode:

pip install -e ".[dev]"
python examples/quick_start.py
python examples/manual_split.py
python examples/deconvolution_demo.py
python examples/full_workflow.py

For headless servers (no display) set the matplotlib backend:

MPLBACKEND=Agg python examples/full_workflow.py