開発版ドキュメント · 0.2.7 12b76261 · 入門例の検証対象 0.2.4 · 版情報 · 既知の制限

イベントカタログと時系列データを対応付ける

マルチメッセンジャー解析や検出器の特性評価ワークフローでは、離散イベントカタログ（突発性グリッチ、環境トリガー、補助環境トリガーなど）と多チャンネル連続時系列を頻繁に対応付けます。カタログのタイムスタンプが厳密であると仮定するのではなく、高レートの環境チャンネルを用いてイベント時刻を精密化し、異種センサを高精度に整列させます。

このチュートリアルで達成できること：

1. 注入された遅延、データ欠損、および端点イベントを含む、8イベントのカタログと3チャンネルの時系列（環境、センサ、変位）を生成します。

2. gwpy.segments.Segment および SegmentTable を使用して探索窓を管理します。

3. 正の環境パルス上で候補タイムスタンプを精密化し、負のトリガーや非ピークトリガーを除外します。

4. 半開区間によるチャンネル固有窓（[start, end)）を抽出し、異種物理単位を混在させることなく有効なイベントを整列します。

5. 構造化カタログ、整列された .npz 配列を出力し、すべてのカタログ完全性チェックを検証します。

データ種別: 合成多チャンネル時系列（600秒、128 Hz）およびイベントカタログ。

環境のセットアップ

import json
import os
import platform
import tempfile
from pathlib import Path

from IPython.display import display
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from astropy import units as u
from scipy import signal

import gwexpy
from gwexpy.timeseries import TimeSeries, TimeSeriesDict
from gwpy.segments import Segment
from gwpy.table import Table

# Define output directory
output_dir_env = os.environ.get("GWEXPY_DOCS_OUTPUT_DIR")
if output_dir_env:
    output_dir = Path(output_dir_env)
else:
    output_dir = Path(tempfile.mkdtemp(prefix="gwexpy-t2-"))
output_dir.mkdir(parents=True, exist_ok=True)
(output_dir / "tables").mkdir(exist_ok=True)
(output_dir / "figures").mkdir(exist_ok=True)
print(f"Artifacts will be written to: {output_dir}")

Artifacts will be written to: /tmp/gwexpy-t2-4bhugutn

カタログと複数チャンネルのデータ契約

duration_s = 600.0
fs = 128.0
dt_s = 1.0 / fs
n_samples = int(duration_s * fs)
gps_t0 = 1400000000.0
time_s = np.arange(n_samples) * dt_s

# Catalog: 8 events E00-E07
# Nominal target offsets: [1, 60, 120, 180, 240, 360, 480, 599] s
target_offsets = np.array([1.0, 60.0, 120.0, 180.0, 240.0, 360.0, 480.0, 599.0])
sample_shifts = np.array([2, -3, 1, 0, -2, 3, -1, 2])
truth_gps = gps_t0 + target_offsets
truth_samples = (target_offsets * fs).astype(int)
catalog_gps = gps_t0 + target_offsets + (sample_shifts * dt_s)
catalog_samples = truth_samples + sample_shifts

catalog_df = pd.DataFrame({
    "event_id": [f"E{i:02d}" for i in range(8)],
    "catalog_gps_s": catalog_gps,
    "catalog_sample": catalog_samples,
    "truth_gps_s": truth_gps,
    "truth_sample": truth_samples,
    "target_offset_s": target_offsets,
})
catalog_path = output_dir / "tables/catalog.csv"
catalog_df.to_csv(catalog_path, index=False)
print(f"Saved initial catalog to {catalog_path}")

# Generate waveforms: witness [V], sensor [V], displacement [m]
rng = np.random.default_rng(2026091602)
val_wit = rng.normal(0, 0.001, n_samples)
val_sen = rng.normal(0, 0.0005, n_samples)
val_disp = rng.normal(0, 1e-12, n_samples)

pulse_width_s = 0.03
pulse_samples = int(pulse_width_s * fs)

# Injected lags: sensor +4 samples, displacement +8 samples
sensor_lag_samples = 4
disp_lag_samples = 8

for i, t_off in enumerate(target_offsets):
    # E05 is a negative control: catalog entry with NO physical pulse in witness
    if i == 5:
        continue
    c_idx = int(t_off * fs)
    # Gaussian pulse
    rel_t = np.arange(-3 * pulse_samples, 3 * pulse_samples + 1)
    g_pulse = np.exp(-0.5 * (rel_t / (pulse_samples / 2.0))**2)
    
    # Witness
    w_idx = c_idx + rel_t
    valid_w = (w_idx >= 0) & (w_idx < n_samples)
    val_wit[w_idx[valid_w]] += 1.0 * g_pulse[valid_w]
    
    # Sensor (lag 4)
    s_idx = c_idx + sensor_lag_samples + rel_t
    valid_s = (s_idx >= 0) & (s_idx < n_samples)
    val_sen[s_idx[valid_s]] += 0.5 * g_pulse[valid_s]
    
    # Displacement (lag 8)
    d_idx = c_idx + disp_lag_samples + rel_t
    valid_d = (d_idx >= 0) & (d_idx < n_samples)
    val_disp[d_idx[valid_d]] += 1e-9 * g_pulse[valid_d]

# E03 sensor gap: inject NaNs into sensor around E03 (178s to 182s)
gap_mask = (time_s >= 178.0) & (time_s <= 182.0)
val_sen[gap_mask] = np.nan

ts_dict = TimeSeriesDict({
    "WITNESS": TimeSeries(val_wit, t0=gps_t0, dt=dt_s * u.s, unit=u.V, channel="WITNESS"),
    "SENSOR": TimeSeries(val_sen, t0=gps_t0, dt=dt_s * u.s, unit=u.V, channel="SENSOR"),
    "DISPLACEMENT": TimeSeries(val_disp, t0=gps_t0, dt=dt_s * u.s, unit=u.m, channel="DISPLACEMENT"),
})
print("Synthesized 3 channels with delays, gaps on E03 sensor, and edge events E00/E07.")

Saved initial catalog to /tmp/gwexpy-t2-4bhugutn/tables/catalog.csv
Synthesized 3 channels with delays, gaps on E03 sensor, and edge events E00/E07.

ウィットネスパルスを用いた時刻の精密化

refined_records = []

for _, row in catalog_df.iterrows():
    eid = row["event_id"]
    c_gps = row["catalog_gps_s"]
    c_sample = int(row["catalog_sample"])
    t_gps = row["truth_gps_s"]
    t_sample = int(row["truth_sample"])
    
    # Search window: catalog +- 0.25 s
    w_start = max(0, c_sample - int(0.25 * fs))
    w_end = min(n_samples, c_sample + int(0.25 * fs))
    sub_wit = val_wit[w_start:w_end]
    
    # Find positive peaks (height >= 0.5V, prominence >= 0.25V)
    peaks, props = signal.find_peaks(sub_wit, height=0.5, prominence=0.25, distance=int(0.1 * fs))
    
    if len(peaks) > 0:
        # If multiple, take largest peak (or earliest if equal)
        best_p = peaks[np.argmax(props["peak_heights"])]
        ref_sample = w_start + best_p
        ref_gps = gps_t0 + ref_sample * dt_s
        status = "refined"
    else:
        ref_sample = c_sample
        ref_gps = c_gps
        status = "no_peak"
        
    delta_correction = (ref_sample - c_sample) * dt_s
    delta_error = (ref_sample - t_sample) * dt_s
        
    refined_records.append({
        "event_id": eid,
        "catalog_sample": c_sample,
        "refined_sample": ref_sample,
        "truth_sample": t_sample,
        "catalog_gps_s": c_gps,
        "refined_gps_s": ref_gps,
        "truth_gps_s": t_gps,
        "delta_correction_s": delta_correction,
        "delta_error_s": delta_error,
        "refine_status": status,
    })

refined_df = pd.DataFrame(refined_records)
refined_path = output_dir / "tables/refined_events.csv"
refined_df.to_csv(refined_path, index=False)
print("Refined event timestamps:")
display(refined_df)

Refined event timestamps:

	event_id	catalog_sample	refined_sample	truth_sample	catalog_gps_s	refined_gps_s	truth_gps_s	delta_correction_s	delta_error_s	refine_status
0	E00	130	128	128	1.400000e+09	1.400000e+09	1.400000e+09	-0.015625	0.000000	refined
1	E01	7677	7680	7680	1.400000e+09	1.400000e+09	1.400000e+09	0.023438	0.000000	refined
2	E02	15361	15360	15360	1.400000e+09	1.400000e+09	1.400000e+09	-0.007812	0.000000	refined
3	E03	23040	23040	23040	1.400000e+09	1.400000e+09	1.400000e+09	0.000000	0.000000	refined
4	E04	30718	30720	30720	1.400000e+09	1.400000e+09	1.400000e+09	0.015625	0.000000	refined
5	E05	46083	46083	46080	1.400000e+09	1.400000e+09	1.400000e+09	0.000000	0.023438	no_peak
6	E06	61439	61440	61440	1.400000e+09	1.400000e+09	1.400000e+09	0.007812	0.000000	refined
7	E07	76674	76672	76672	1.400001e+09	1.400001e+09	1.400001e+09	-0.015625	0.000000	refined

カタログ時刻と精密化時刻の可視化

fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(time_s[:int(65 * fs)], val_wit[:int(65 * fs)], label="Witness Waveform", color="navy", lw=0.8)
for _, r in refined_df.head(2).iterrows():
    ax.axvline(r["catalog_gps_s"] - gps_t0, color="gray", ls="--", label=f"{r['event_id']} Catalog" if r['event_id'] == 'E00' else None)
    ax.axvline(r["refined_gps_s"] - gps_t0, color="crimson", ls="-", label=f"{r['event_id']} Refined" if r['event_id'] == 'E00' else None)

ax.set_xlim(0, 65)
ax.set_xlabel("Time from Start [s]")
ax.set_ylabel("Witness [V]")
ax.set_title("Catalog vs Refined Peak Timestamps (First 65 seconds)")
ax.grid(True, alpha=0.3)
ax.legend(loc="upper right")
fig.tight_layout()

fig_path_cat = output_dir / "figures/catalog_refinement.png"
fig.savefig(fig_path_cat, dpi=120)
display(fig)
plt.close(fig)
print(f"Saved catalog refinement plot to {fig_path_cat}")

Saved catalog refinement plot to /tmp/gwexpy-t2-4bhugutn/figures/catalog_refinement.png

チャンネル別半開区間ウィンドウの切り出し

# Windows config: witness [-0.25, 0.5), sensor [-1.0, 2.0), displacement [-2.0, 3.0)
win_configs = {
    "WITNESS": (-0.25, 0.50),
    "SENSOR": (-1.00, 2.00),
    "DISPLACEMENT": (-2.00, 3.00),
}

channel_records = []
aligned_data = {"WITNESS": [], "SENSOR": [], "DISPLACEMENT": []}
valid_event_ids = {"WITNESS": [], "SENSOR": [], "DISPLACEMENT": []}

for _, r in refined_df.iterrows():
    eid = r["event_id"]
    ref_s = r["refined_sample"]
    ref_status = r["refine_status"]
    
    for ch_name, (pre_s, post_s) in win_configs.items():
        ts_ch = ts_dict[ch_name]
        req_start = ref_s + int(pre_s * fs)
        req_end = ref_s + int(post_s * fs)
        req_n = req_end - req_start
        
        act_start = max(0, req_start)
        act_end = min(n_samples, req_end)
        act_n = act_end - act_start
        
        # Check boundary & gap conditions
        if ref_status == "no_peak":
            status = "no_peak"
            seg_val = None
        elif act_n < req_n:
            status = "partial"
            seg_val = None
        else:
            seg_val = ts_ch.value[act_start:act_end]
            if not np.all(np.isfinite(seg_val)):
                status = "gap"
                seg_val = None
            else:
                status = "complete"
                aligned_data[ch_name].append(seg_val)
                valid_event_ids[ch_name].append(eid)
                
        peak_val = float(np.nanmax(np.abs(seg_val))) if seg_val is not None else np.nan
        rms_val = float(np.sqrt(np.nanmean(seg_val**2))) if seg_val is not None else np.nan
        
        channel_records.append({
            "event_id": eid,
            "channel": ch_name,
            "requested_samples": req_n,
            "actual_samples": act_n,
            "peak_val": peak_val,
            "rms_val": rms_val,
            "unit": str(ts_ch.unit),
            "status": status,
        })

ch_df = pd.DataFrame(channel_records)
ch_path = output_dir / "tables/event_channels.csv"
ch_df.to_csv(ch_path, index=False)
print(f"Channel window summary ({len(ch_df)} rows = 8 events x 3 channels):")
display(ch_df.head(12))

# Save aligned arrays to .npz (allow_pickle=False)
for ch_name, (pre_s, post_s) in win_configs.items():
    rel_time = np.arange(int((post_s - pre_s) * fs)) * dt_s + pre_s
    arr = np.array(aligned_data[ch_name])
    np.savez(
        output_dir / f"tables/aligned_{ch_name.lower()}.npz",
        relative_time_s=rel_time,
        waveforms=arr,
        event_ids=np.array(valid_event_ids[ch_name]),
        sample_rate_hz=fs,
    )

Channel window summary (24 rows = 8 events x 3 channels):

	event_id	channel	requested_samples	actual_samples	peak_val	rms_val	unit	status
0	E00	WITNESS	96	96	9.991960e-01	1.664244e-01	V	complete
1	E00	SENSOR	384	384	4.989713e-01	4.160944e-02	V	complete
2	E00	DISPLACEMENT	640	512	NaN	NaN	m	partial
3	E01	WITNESS	96	96	9.994821e-01	1.665012e-01	V	complete
4	E01	SENSOR	384	384	4.995849e-01	4.159562e-02	V	complete
5	E01	DISPLACEMENT	640	640	9.978954e-10	6.440016e-11	m	complete
6	E02	WITNESS	96	96	1.000068e+00	1.665091e-01	V	complete
7	E02	SENSOR	384	384	4.999338e-01	4.161557e-02	V	complete
8	E02	DISPLACEMENT	640	640	9.989395e-10	6.440815e-11	m	complete
9	E03	WITNESS	96	96	9.989518e-01	1.662120e-01	V	complete
10	E03	SENSOR	384	384	NaN	NaN	V	gap
11	E03	DISPLACEMENT	640	640	1.000935e-09	6.448848e-11	m	complete

整列波形と平均応答の可視化

fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=False)

for ax, ch_name in zip(axes, ["WITNESS", "SENSOR", "DISPLACEMENT"]):
    pre_s, post_s = win_configs[ch_name]
    rel_time = np.arange(int((post_s - pre_s) * fs)) * dt_s + pre_s
    arr = np.array(aligned_data[ch_name])
    for w in arr:
        ax.plot(rel_time, w, color="gray", alpha=0.5, lw=0.8)
    if len(arr) > 0:
        ax.plot(rel_time, np.mean(arr, axis=0), color="crimson", lw=1.8, label=f"Mean (N={len(arr)})")
    ax.set_ylabel(f"{ch_name} [{ts_dict[ch_name].unit}]")
    ax.set_xlabel("Relative Time to Refined Event [s]")
    ax.grid(True, alpha=0.3)
    ax.legend(loc="upper right")

fig.suptitle("Aligned Multi-Channel Event Waveforms")
fig.tight_layout()

fig_path_align = output_dir / "figures/aligned_waveforms.png"
fig.savefig(fig_path_align, dpi=120)
display(fig)
plt.close(fig)

# Figure 3: Mean response
fig2, ax = plt.subplots(figsize=(8, 4))
rel_t_wit = np.arange(int(0.75 * fs)) * dt_s - 0.25
mean_wit = np.mean(aligned_data["WITNESS"], axis=0)
ax.plot(rel_t_wit, mean_wit, color="darkblue", lw=2, label="Witness Mean Response")
ax.set_xlabel("Relative Time [s]")
ax.set_ylabel("Witness [V]")
ax.set_title("Coherent Witness Mean Response across Complete Events")
ax.grid(True, alpha=0.3)
ax.legend()
fig2.tight_layout()
fig_path_mean = output_dir / "figures/mean_response.png"
fig2.savefig(fig_path_mean, dpi=120)
display(fig2)
plt.close(fig2)
print("Saved aligned waveform and mean response figures.")

Saved aligned waveform and mean response figures.

検証とカタログ整合性検査

# 1. Recovered lag check (sensor lag = 4, disp lag = 8)
# Compute peak delay on complete inner event E01
sub_wit = ts_dict["WITNESS"].value[int(59 * fs):int(61 * fs)]
sub_sen = ts_dict["SENSOR"].value[int(59 * fs):int(61 * fs)]
sub_disp = ts_dict["DISPLACEMENT"].value[int(59 * fs):int(61 * fs)]

p_wit = np.argmax(sub_wit)
p_sen = np.argmax(sub_sen)
p_disp = np.argmax(sub_disp)

obs_sen_lag = float(p_sen - p_wit)
obs_disp_lag = float(p_disp - p_wit)
print(f"Observed lags: sensor={obs_sen_lag} samples (expected 4), disp={obs_disp_lag} samples (expected 8)")

# 2. Invalid catalog negative checks
# Check duplicate ID rejection
def validate_catalog(df):
    if len(df["event_id"].unique()) != len(df):
        raise ValueError("Duplicate event IDs detected")
    if df["catalog_gps_s"].isna().any():
        raise ValueError("NaN timestamps detected in catalog")
    return True

assert validate_catalog(catalog_df) is True
try:
    bad_cat = pd.DataFrame({"event_id": ["E01", "E01"], "catalog_gps_s": [1.0, 2.0]})
    validate_catalog(bad_cat)
    cat_check_passed = False
except ValueError:
    cat_check_passed = True
assert cat_check_passed is True
print("Catalog negative checks passed.")

# 3. 8x3 correspondence, ID set uniqueness, units and sample rates
assert len(ch_df) == 24, "Must have exactly 24 event-channel rows (8 events x 3 channels)"
assert set(ch_df["event_id"].unique()) == set(catalog_df["event_id"].unique()), "Event ID set must match catalog"
assert len(catalog_df["event_id"].unique()) == 8, "Must have 8 unique event IDs"
assert (ch_df.groupby("event_id").size() == 3).all(), "Every event must map to exactly 3 channels"
channel_units = {ch: str(ts_dict[ch].unit) for ch in ["WITNESS", "SENSOR", "DISPLACEMENT"]}
channel_rates = {ch: float(ts_dict[ch].sample_rate.value) for ch in ["WITNESS", "SENSOR", "DISPLACEMENT"]}
assert channel_units == {"WITNESS": "V", "SENSOR": "V", "DISPLACEMENT": "m"}
assert all(r == fs for r in channel_rates.values())
print(f"Verified 8x3 correspondence, ID uniqueness, units ({channel_units}), and sample rates ({channel_rates}).")

Observed lags: sensor=4.0 samples (expected 4), disp=8.0 samples (expected 8)
Catalog negative checks passed.
Verified 8x3 correspondence, ID uniqueness, units ({'WITNESS': 'V', 'SENSOR': 'V', 'DISPLACEMENT': 'm'}), and sample rates ({'WITNESS': 128.0, 'SENSOR': 128.0, 'DISPLACEMENT': 128.0}).

解析設定と検証メトリクスの保存

settings = {
    "tutorial_id": "T2",
    "data_kind": "synthetic",
    "seed": 2026091602,
    "gps_t0_s": gps_t0,
    "sample_rate_hz": fs,
    "duration_s": duration_s,
    "channel_units": {"WITNESS": "V", "SENSOR": "V", "DISPLACEMENT": "m"},
    "analysis_parameters": {"search_window_s": 0.25, "sensor_lag": 4, "disp_lag": 8},
    "python_version": platform.python_version(),
}
with open(output_dir / "analysis-settings.json", "w", encoding="utf-8") as f:
    json.dump(settings, f, indent=2)

metrics = {
    "status": "passed",
    "data_kind": "synthetic",
    "checks": {
        "event_id_preservation": {
            "observed": {
                "total_rows": int(len(ch_df)),
                "unique_events": int(len(ch_df["event_id"].unique())),
                "unique_pairs_count": int(len(set(zip(ch_df["event_id"], ch_df["channel"])))),
                "has_duplicate_pairs": bool(ch_df.duplicated(subset=["event_id", "channel"]).any()),
                "matches_cartesian_product": bool(
                    set(zip(ch_df["event_id"], ch_df["channel"])) ==
                    {(eid, ch) for eid in catalog_df["event_id"].unique() for ch in ["WITNESS", "SENSOR", "DISPLACEMENT"]}
                ),
            },
            "criterion": "exactly 24 unique (event_id, channel) pairs matching Cartesian product of 8 catalog events x 3 channels",
            "passed": bool(
                len(ch_df) == 24 and
                not ch_df.duplicated(subset=["event_id", "channel"]).any() and
                set(zip(ch_df["event_id"], ch_df["channel"])) ==
                {(eid, ch) for eid in catalog_df["event_id"].unique() for ch in ["WITNESS", "SENSOR", "DISPLACEMENT"]}
            ),
        },
        "event_refinement_accuracy": {
            "observed": float(np.max(np.abs(refined_df[refined_df["refine_status"] == "refined"]["delta_error_s"]))),
            "criterion": "truth estimation error <= 1 / fs (0.0078125 s)",
            "passed": bool(np.max(np.abs(refined_df[refined_df["refine_status"] == "refined"]["delta_error_s"])) <= (1.0 / fs)),
            "max_correction_offset_s": float(np.max(np.abs(refined_df[refined_df["refine_status"] == "refined"]["delta_correction_s"]))),
        },
        "event_lag_recovery": {
            "observed": [obs_sen_lag, obs_disp_lag],
            "criterion": "sensor lag == 4 and disp lag == 8 within 1 sample",
            "passed": bool(abs(obs_sen_lag - 4.0) <= 1.0 and abs(obs_disp_lag - 8.0) <= 1.0),
        },
        "event_window_half_open": {
            "observed": int(ch_df[ch_df["status"] == "complete"]["requested_samples"].iloc[0]),
            "criterion": "half-open window length matches (post - pre) * fs",
            "passed": bool(((ch_df["requested_samples"] == ch_df["actual_samples"]) | (ch_df["status"] != "complete")).all()),
        },
        "event_edge_and_gap_status": {
            "observed": list(ch_df[ch_df["status"].isin(["partial", "gap"])]["event_id"].unique()),
            "criterion": "partial edge events (E00, E07) and gap (E03) flagged",
            "passed": bool("E03" in list(ch_df[ch_df["status"] == "gap"]["event_id"]) and "E00" in list(ch_df[ch_df["status"] == "partial"]["event_id"])),
        },
        "event_alignment_no_unit_mixing": {
            "observed": {
                "units": channel_units,
                "sample_rates_hz": channel_rates,
            },
            "criterion": "channels aligned independently with strict units (WITNESS: V, SENSOR: V, DISPLACEMENT: m) and 128 Hz sampling",
            "passed": bool(
                ts_dict["WITNESS"].unit == u.V and
                ts_dict["SENSOR"].unit == u.V and
                ts_dict["DISPLACEMENT"].unit == u.m and
                all(ts_dict[ch].sample_rate.value == fs for ch in ["WITNESS", "SENSOR", "DISPLACEMENT"]) and
                len(aligned_data["WITNESS"]) > 0 and len(aligned_data["DISPLACEMENT"]) > 0
            ),
        },
        "event_catalog_invalid": {
            "observed": cat_check_passed,
            "criterion": "duplicate or invalid catalog entries rejected",
            "passed": bool(cat_check_passed is True),
        },
    },
}
with open(output_dir / "validation-metrics.json", "w", encoding="utf-8") as f:
    json.dump(metrics, f, indent=2)

print("Settings and metrics saved successfully.")

Settings and metrics saved successfully.