バックトレース (`data.trace`) — 実験的¶

data.trace は、6D 位相空間グリッドから粒子を作り、到達確率、backward trace、forward trace をまとめて扱う workflow API です。新しい解析コードでは data.trace.forward() / data.trace.backward() / data.trace.both() を使ってください。

従来の data.backtrace は、単一粒子や任意の Particle 配列を直接渡したい既存コード向けの低水準 API として残っています。このページでは本線を data.trace にし、data.backtrace の経路は折りたたみ内にまとめています。

要件: バックトレースは外部の vdist-solver-fortran パッケージ (vdsolverf) に依存します。pip install vdist-solver-fortran でインストールしてください。未インストールの環境では data.trace.* の呼び出し時に ImportError が出ます。

いつ使うか¶

ある観測点に到達する粒子の 位相空間分布 を見たい
到達粒子を backward / forward の軌跡 として追いたい
到達粒子の エネルギースペクトル を描きたい
確率と軌跡を同じ粒子集合から計算して、確率を alpha に使いたい

バックトレースは保存済みフィールドを使って ODE を解くため、大きな max_step や細かい位相空間グリッドでは時間がかかります。HPC ノードに処理を任せたい場合は、リモート実行と組み合わせてください（後述）。

単位に注意: data.trace の入力は EMSES シミュレーション単位です。SI 値から指定する場合は data.unit で変換してから渡します。

クイックスタート¶

import emout

data = emout.Emout("output_dir")

vx_scan = (data.unit.v.trans(-3e5), data.unit.v.trans(3e5), 64)
vz_scan = (data.unit.v.trans(-3e5), data.unit.v.trans(3e5), 64)

trace = data.trace.forward(
    x=20.0, y=32.0, z=40.0,
    vx=vx_scan,
    vy=0.0,
    vz=vz_scan,
    ispec=0,
    get_trace=True,
    get_probabilities=True,
    max_step=10000,
    n_threads=8,
)

trace.plot("vx", "vz", cmap="viridis")       # 到達確率 heatmap
trace.plot_traces("x", "z")                  # 確率 alpha 付き軌跡
trace.probabilities.plot_energy_spectrum(scale="log")

スカラー値だけを渡すと、位相空間グリッドのサイズ 1 軸として扱われます。確率計算なしで 1 本の軌跡だけ欲しい場合も data.trace で書けます:

single = data.trace.backward(
    x=20.0, y=32.0, z=40.0,
    vx=data.unit.v.trans(1.0e5),
    vy=0.0,
    vz=data.unit.v.trans(-2.0e5),
    ispec=0,
    get_trace=True,
    get_probabilities=False,
)

single.plot_traces("t", "x")
single.traces.xvz.plot()

入力単位の約束¶

data.trace に渡す x / y / z / vx / vy / vz、dt、probability_dt は すべて EMSES シミュレーション単位です。emout は SI から自動変換しません。

SI 値から指定したい場合は、呼び出し前に data.unit で EMSES 単位へ変換してください:

position = (
    data.unit.length.trans(0.20),  # m -> EMSES length
    data.unit.length.trans(0.32),
    data.unit.length.trans(0.40),
)
vx_scan = (
    data.unit.v.trans(-3.0e5),     # m/s -> EMSES velocity
    data.unit.v.trans(3.0e5),
    64,
)
vz_scan = (
    data.unit.v.trans(-3.0e5),
    data.unit.v.trans(3.0e5),
    64,
)

data.unit は plasma.inp の !!key dx=...,to_c=... ヘッダー（または plasma.toml の [meta.unit_conversion]）がある場合にだけ利用できます。単位変換 metadata がない場合は、EMSES 単位の値を直接渡してください。

TraceResult.phases、TraceResult.particles、trace.traces.positions_list などの配列も EMSES 単位のまま保持されます。一方で、trace.plot() / trace.plot_traces() / trace.traces.xz.plot() は unit metadata がある場合に既定で SI 単位へ変換して表示します（use_si=False で EMSES 単位表示）。

dt と probability_dt は EMSES 時間単位の非負の刻み幅です。None の場合は abs(data.inp.dt) を使います。data.trace.backward() は solver に +dt、data.trace.forward() は内部で符号を反転して -dt を渡します。data.trace.both() は同じ dt から backward に +dt、forward に -dt を使います。到達確率計算は probability_dt を backward 側の符号で使います。負の値を渡すと ValueError になります。

電場・磁場を無効化する¶

use_electric_field と use_magnetic_field はキーワード専用引数で、どちらも既定値は True です。False を指定すると、対応する場の読み込みと粒子への作用を無効化します。False を指定するには vdist-solver-fortran >= 1.7.0 が必要です。

引数	`False` で無効になる場
`use_electric_field`	空間電荷と蓄積電荷による電場
`use_magnetic_field`	出力された磁場と `wc` から作る背景磁場

例えば、電場だけを使って確率と両方向の軌跡を計算するには:

trace = data.trace.both(
    x=20.0, y=32.0, z=40.0,
    vx=vx_scan, vy=0.0, vz=vz_scan,
    get_trace=True,
    get_probabilities=True,
    use_magnetic_field=False,
)

両方を False にすることもできます。設定は確率計算とすべての軌跡計算に共通で、境界・衝突の扱いは通常どおりです。

data.trace.forward() / backward() / both() と、低水準 data.backtrace のすべての get_* メソッドで使えます。Dask リモート実行や parallel="mpi" / parallel="srun" でも同じ指定が渡されます。省略時は従来どおり両方の場を使い、旧版 solver にも対応します。

Workflow API: `data.trace`¶

data.trace.backward() / data.trace.forward() / data.trace.both() は常に TraceResult を返します。要求しなかった payload は None です。

確率と軌跡をまとめて得る¶

trace = data.trace.forward(
    x=20.0, y=32.0, z=40.0,
    vx=vx_scan,
    vy=0.0,
    vz=vz_scan,
    get_trace=True,
    get_probabilities=True,
)

trace.probabilities        # ProbabilityResult
trace.forward_traces       # MultiBacktraceResult
trace.backward_traces      # None
trace.alpha                # np.clip(trace.probabilities.probabilities, 0, 1)

trace.plot("vx", "vz")     # 到達確率 heatmap
trace.plot_traces("x", "z")

trace.plot() は確率がある場合は trace.probabilities.pair(...).plot()、確率がない場合は trace.plot_traces() に自動で振り分けます。明示したい場合は kind="probability" / kind="trace" を渡します。

軌跡だけを得る¶

get_probabilities=False にすると、確率計算を省略して phase-space grid から粒子だけを作り、軌跡だけを得ます。この場合 trace.probabilities と trace.alpha は None で、plot_traces() は一様な alpha で描画します。

trace = data.trace.forward(
    x=20.0, y=32.0, z=40.0,
    vx=vx_scan,
    vy=0.0,
    vz=vz_scan,
    get_trace=True,
    get_probabilities=False,
)

trace.forward_traces.xz.plot(alpha=0.3)
trace.plot_traces("x", "z", alpha=0.3)

backward と forward を同時に得る¶

both() は同じ phase-space grid から backward / forward の両方を計算します。確率も要求した場合、確率計算は 1 回だけです。

trace = data.trace.both(..., get_trace=True)
trace.plot_traces("x", "z")
trace.plot_traces("x", "z", direction="backward")
trace.plot_traces("x", "z", direction="forward")

both() の結果では、plot_traces() の direction を省略すると backward / forward の両方を同じ axes に重ねます。片方向だけを描く場合は direction="backward" または direction="forward" を指定してください。trace.traces は単一の trajectory payload を返す属性なので、both() の結果から直接参照した場合は引き続き曖昧性エラーになります。

3D に重ねる¶

plot3d() は PyVista plotter を返します。both() の結果では direction を省略すると両方向を描きます。既存の 3D field / boundary 表示に重ねたい場合は、その plotter を渡します。

plotter = data.phisp[-1].plot3d(mode="slice", show=False)
trace.plot3d(plotter=plotter, tube_radius=0.05, show=True)

小さい確率の軌跡を見えるようにするには min_alpha=0.08、相対的な分布を表示するには alpha_scale="normalize" / "log" を指定します。direction="both"、方向別の色、表示単位での clip_bounds、目盛とカメラの指定も利用できます。境界の電位着色と組み合わせた保存例は PyVista 可視化を参照してください。

結果の描画¶

2D ヒートマップへの射影¶

trace.probabilities.pair(var1, var2) は選ばれていない 4 軸を台形則で積分し、HeatmapData を返します。trace.plot(var1, var2) はこの操作を短く書く入口です。

trace.plot("vx", "vz", cmap="viridis")
trace.probabilities.xvx.plot()
trace.probabilities.yz.plot(cmap="plasma")

HeatmapData.plot() はカラーバー付きの pcolormesh を描き、軸ラベルには SI 単位が付きます（use_si=False でグリッド単位に戻せます）。追加のキーワード引数はすべて pcolormesh にそのまま渡されるので、vmin / vmax や norm=LogNorm(...) で配色を制御できます。

エネルギースペクトル¶

エネルギースペクトルは確率 payload 側の ProbabilityResult から描きます。

trace.probabilities.plot_energy_spectrum(scale="log", energy_bins=80)
hist, bin_edges = trace.probabilities.energy_spectrum(energy_bins=80)

内部で plasma.inp の wp や光電流設定から基準数密度 n0 を計算し、各位相点の確率に掛けてから積分しています（nflag_emit == 2 の光電子は path と curf から導出）。

リモート実行との連携¶

data.trace は Emout ファサード経由で remote_open_kwargs を持ち、emout server が動いていれば既定でワーカー上で実行されます（結果は RemoteTraceResult プロキシとして返ります）。ワーカー上にキャッシュされているので、プロットパラメータを変えても再計算は走りません。

from emout.distributed import remote_figure

trace = data.trace.forward(
    x=20.0, y=32.0, z=40.0,
    vx=vx_scan,
    vy=0.0,
    vz=vz_scan,
    get_trace=True,
)

with remote_figure():
    trace.plot("vx", "vz", cmap="viridis")

with remote_figure():
    trace.plot_traces("x", "z")

trace.drop()   # 不要になったらワーカーメモリを解放

明示的な remote 記法を使いたいときは data.remote().trace... に切り替えます:

from emout.distributed import remote_scope, remote_figure

with remote_scope():
    rdata = data.remote()
    trace = rdata.trace.forward(
        x=20.0, y=32.0, z=40.0,
        vx=vx_scan,
        vy=0.0,
        vz=vz_scan,
        get_trace=True,
    )

    with remote_figure():
        trace.plot("vx", "vz")
        trace.plot_traces("x", "z")

リモート実行全体の仕組み・環境変数・サーバー起動についてはリモート実行ガイドを参照してください。

`fetch()` でローカルに取り出す¶

matplotlib で細かくカスタマイズしたい（独自アノテーション、共通カラーバー、subplot への差し込み、など）場合は、fetch() で結果を小さなローカルオブジェクトに変換できます:

local_trace = trace.fetch()
heatmap = local_trace.probabilities.vxvz
fig, ax = plt.subplots()
heatmap.plot(ax=ax, cmap="plasma")
ax.axhline(y=0, color="red", linestyle="--")

既存コード向け: 低水準 API `data.backtrace` を表示

低水準 API: data.backtrace

data.backtrace は、data.trace が内部で使っている BacktraceWrapper です。1 本の軌跡を明示的な position / velocity で計算したい場合、vdsolverf.core.Particle を直接渡したい場合、または solver に渡す dt の符号を直接制御したい場合に使います。新しい位相空間 workflow は data.trace を優先してください。

単一粒子: get_backtrace

position = (20.0, 32.0, 40.0)
velocity = (
    data.unit.v.trans(1.0e5),
    0.0,
    data.unit.v.trans(-2.0e5),
)

bt = data.backtrace.get_backtrace(position, velocity, ispec=0, max_step=50000)

bt.tx.plot()                 # = bt.pair("t", "x")
bt.xvz.plot()                # = bt.pair("x", "vz")
bt.yz.plot(color="black")    # 軌道の yz 投影

bt.ts、bt.probability、bt.positions、bt.velocities は EMSES 単位の配列です。XYData.plot() はデフォルトで SI 単位に変換し、軸ラベルも自動生成します（use_si=False で EMSES 単位のまま表示できます）。

低水準 API で forward trace を行う場合は dt=-abs(data.inp.dt) を渡します。 position / velocity / dt はいずれも EMSES 単位で、描画時の SI 変換とは独立です。高水準の data.trace.forward() は非負の刻み幅を受け取り、内部で符号を変えます。

複数粒子: get_backtraces

import numpy as np

positions = np.array([[20, 32, 40], [21, 32, 40], [22, 32, 40]], dtype=float)
velocities = np.zeros_like(positions)
velocities[:, 0] = data.unit.v.trans(1.0e5)

many = data.backtrace.get_backtraces(positions, velocities, ispec=0)
many.xz.plot(alpha=0.5)
many.sample(50, random_state=0).tvx.plot()

positions と velocities は対応する粒子ごとの (N, 3) 配列です。Cartesian product の位相空間グリッドを作りたい場合は data.trace を使ってください。

Particle オブジェクトを直接渡す

from vdsolverf.core import Particle

particles = [Particle(p, v) for p, v in zip(positions, velocities)]
many = data.backtrace.get_backtraces_from_particles(particles, ispec=0)

ProbabilityResult から得られた粒子リストをそのまま渡して、確率が高い粒子だけ軌跡を描く、といった連結もできます:

result = data.backtrace.get_probabilities(...)
bt = data.backtrace.get_backtraces_from_particles(result.particles, ispec=0)
bt.xz.plot(alpha=np.clip(result.probabilities, 0, 1))

到達確率: get_probabilities

data.trace は内部で data.backtrace.get_probabilities(...) を呼び、ProbabilityResult を trace.probabilities として保持します。低水準 API を直接使う場合は次の形です:

result = data.backtrace.get_probabilities(
    x=20.0, y=32.0, z=40.0,
    vx=vx_scan,
    vy=0.0,
    vz=vz_scan,
    ispec=0,
    max_step=10000,
    n_threads=8,
)

result.vxvz.plot(cmap="viridis")
result.plot_energy_spectrum(scale="log")

MPI backend

parallel="mpi" / parallel="srun" は低水準 get_probabilities() の backend オプションです。data.trace にも **kwargs として渡せます。

trace = data.trace.forward(
    x=20.0, y=32.0, z=40.0,
    vx=vx_scan,
    vy=0.0,
    vz=vz_scan,
    max_step=10000,
    parallel="srun",
    ntasks=8,
    n_threads=2,
    cpus_per_task=2,
)

HPC の backend 診断¶

計算を実行するノード・worker 上で data.trace.backend_info() または data.backtrace.backend_info() を呼ぶと、解決されたライブラリパス、存在、backend バージョン、 CPU feature、OMP_NUM_THREADS を確認できます。既定ではライブラリをロードしません。

info = data.trace.backend_info()
print(info["library_path"], info["exists"])
probe = data.trace.backend_info(library_path="/path/to/backend.so", check_load=True, timeout=10)
print(probe["loadable"], probe["load_error"])

check_load=True は別の Python プロセスでロードだけを試し、失敗理由・終了シグナル・タイムアウトを返します。solver は実行しません。ロード成功は、計算中の CPU 命令や OpenMP の互換性を保証しません。SIGILL は ISA / ビルド対象、libifport.so 等の不足は対応するコンパイラ runtime を確認する手掛かりです。-xHost 等で別 CPU 向けに作ったライブラリは対象ノード向けの共通命令セットで再ビルドし、OpenMP の異常は serial ビルドと最小ケースで切り分けてください。診断は呼出し元で動くため、リモート実行時は実際の worker 環境で確認します。

検証したライブラリを trace に使うときも、同じパスを明示します。

trace = data.trace.both(
    x=20.0, y=32.0, z=40.0, vx=0.0, vy=0.0, vz=-1.0,
    get_trace=True, library_path="/path/to/backend.so",
)

バックトレース (data.trace) — 実験的¶