粒子の移動と近接関係の更新を別々に測定する
粒子シミュレーションの位置記録から、粒子集合の移動と内部の関係更新を測定します。移動には重心速度、関係更新には近接グラフの構造振幅を使用します。近接グラフとは、一定の距離より近い粒子対をエッジで結んだグラフです。構造振幅は、グラフの時間的な更新量を表します[1]。
次の表は、この記事で検証する状態と観測量の対応です。
| 粒子集合の状態 | 重心速度 | 構造振幅 |
|---|---|---|
| 近接関係を保った平行移動 | 非ゼロ | ゼロ |
| 重心を保った近接関係の更新 | ゼロ | 非ゼロ |
位置の更新と観測量の計算を別々の関数に実装すると、同じ観測処理を異なる生成規則に適用できます。以下では既知の座標系列で観測処理を検証した後、二粒子の相互作用による位置更新にも適用します。
実行環境
検証に使用した環境は次のとおりです。
| 項目 | 値 |
|---|---|
| OS | Debian GNU/Linux 13 (trixie) |
| アーキテクチャ | x86_64 |
| Python | 3.13.5 |
| NumPy | 2.2.4 |
| 座標配列の型 | numpy.float64 |
二次元のユークリッド空間で、2 または 4 粒子の固定集合を扱います。周期境界や乱数を使用せず、各粒子の質量は等しいものとします。
Python 3.13 を使用して仮想環境を作成し、NumPy のバージョンを固定します。
python3.13 -m venv .venv
source .venv/bin/activate
python -m pip install numpy==2.2.4
次のコマンドで使用する処理系とライブラリのバージョンを確認します。
python --version
python -c "import numpy as np; print(np.__version__)"
以下の実行結果は、上記のバージョンが導入済みの Debian 環境で取得したものです。仮想環境を新規作成するコマンドは再現用の手順です。
粒子の位置から近接グラフを構築する
粒子 ID と座標配列
粒子数を $N$ とし、位置を形状 (N, 2) の NumPy 配列で表します。配列の各行は固定した粒子 ID に対応します。観測関数には、空でなく、すべての成分が有限値である浮動小数点配列を渡します。更新前後の配列は同じ形状と粒子 ID の対応を持ち、観測間隔と近接距離には有限の正値を使用します。
4 粒子の初期配置を次のように定めます。
X_0 =
\begin{pmatrix}
-0.5 & -1.5 \\
0.5 & -1.5 \\
-0.5 & 1.5 \\
0.5 & 1.5
\end{pmatrix}
| 行番号 | 粒子 ID | 座標 |
|---|---|---|
| 0 | A | (-0.5, -1.5) |
| 1 | B | (0.5, -1.5) |
| 2 | C | (-0.5, 1.5) |
| 3 | D | (0.5, 1.5) |
時刻間で同じ行番号が同じ粒子を指すようにします。粒子の並べ替えによって ID の対応が失われると、座標の集合が同じでもエッジの追加や削除として計上される場合があります。
距離閾値によるエッジの抽出
離散的な観測時刻を $t$、粒子 $i$ の位置を $\mathbf{x}_{i,t}$、観測上の近接距離を $r_c$ とします。無向エッジの集合を次のように定義します。
E_t = \left\{\{i,j\}\mid i<j,\quad
\|\mathbf{x}_{j,t}-\mathbf{x}_{i,t}\|_2<r_c\right\}
proximity_edges() は、各粒子対の座標差からユークリッド距離を計算します。一次元配列に対して np.linalg.norm() を既定の引数で呼び出すと、2 ノルムが得られます[2]。
近接距離は 2.0 とします。距離が 2.0 未満の粒子対を (i, j) のタプルとして Python の集合に格納します。$i<j$ とすることで、自己ループと逆順の重複を除きます。距離が閾値と等しい粒子対は対象外です。
初期配置では A と B、および C と D の距離がそれぞれ 1.0 であり、それ以外の粒子対の距離は 3.0 以上です。このため、期待するエッジ集合は {(0, 1), (2, 3)} です。
このエッジは観測上の近接関係を表します。相互作用の係数や作用範囲は、粒子の位置を更新する生成規則で別に定めます。
重心速度と構造振幅を計算する
重心の移動速度
粒子の質量が等しい場合、重心は座標の算術平均です。
\mathbf{X}_t = \frac{1}{N}\sum_{i=1}^{N}\mathbf{x}_{i,t}
positions.mean(axis=0) は、各座標成分について粒子方向の平均を計算します。観測間隔を $\Delta t>0$ とすると、連続する二つの位置記録から区間平均の重心速度を求められます。
\mathbf{V}_t = \frac{\mathbf{X}_{t+1}-\mathbf{X}_t}{\Delta t}
centroid_velocity() は、更新前後の位置配列と観測間隔を受け取り、この式を計算します。今回の対象は固定した粒子集合全体です。各時刻で最大連結成分を抽出する処理とは区別します。
エッジ集合の対称差による構造振幅
構造をグラフ $S_t$ として表し、その数値表現を $\phi(S_t)$ とすると、構造振幅は隣接する観測時刻の数値表現間の距離として定義できます[3]。
固定した粒子集合の各無向エッジ候補を一定の順序に並べ、エッジの有無を 0 または 1 で表します。ユークリッド距離を使用すると、構造振幅は次のようになります。
\begin{aligned}
A_{t+1}
&=\|\phi(S_{t+1})-\phi(S_t)\|_2\\
&=\sqrt{|E_{t+1}\triangle E_t|}
\end{aligned}
$\triangle$ は集合の対称差です。どちらか一方だけに含まれる要素を集めます。二値ベクトルの成分差の二乗は、エッジが追加または削除された場合に 1、それ以外の場合に 0 です。このため、ベクトルを明示的に構築せず、対称差の要素数から構造振幅を計算できます。
structural_amplitude() では集合演算 before_edges ^ after_edges によって対称差を求め、その要素数の平方根を返します。
エッジ数の変化量も別に計算します。
\Delta K_{t+1}=|E_{t+1}|-|E_t|
例えば、1 本のエッジが削除され、別の 1 本が追加される場合、エッジ数の変化量は 0、構造振幅は $\sqrt{2}$ です。
観測関数 observe_transition() は、更新前後のエッジ集合、重心速度、対称差の要素数、構造振幅、エッジ数の変化量をまとめて返します。重心速度 $\mathbf{V}_t$ と構造振幅 $A_{t+1}$ は、いずれも時刻 $t$ から $t+1$ までの同一の観測区間に対応します。
既知の座標系列で測定処理を検証する
実験 1：近接関係を維持した平行移動
初期配置の全粒子を右方向へ 1.0 移動させます。観測間隔は 1.0、近接距離は 2.0 とします。
\mathbf{x}_{i,1}=\mathbf{x}_{i,0}+(1,0)
全粒子の位置に同じ変位を加えるため、粒子間の座標差は維持されます。重心は原点から (1, 0) へ移動します。
期待値は次のとおりです。
\mathbf{V}_0=(1,0),\qquad A_1=0,\qquad \Delta K_1=0
更新前後のエッジ集合がともに {(0, 1), (2, 3)} であることも確認します。重心速度が非ゼロでも近接グラフが維持される事例です。
実験 2：重心を維持した近接関係の更新
実験 1 と同じ初期配置から、更新後の位置を次のように設定します。
X_1=\begin{pmatrix}
-1.5 & -0.5 \\
1.5 & -0.5 \\
-1.5 & 0.5 \\
1.5 & 0.5
\end{pmatrix}
近接距離と観測間隔は実験 1 と同じです。更新前後のエッジ集合は次のようになります。
\begin{aligned}
E_0&=\{\{A,B\},\{C,D\}\}\\
E_1&=\{\{A,C\},\{B,D\}\}
\end{aligned}
2 本のエッジが削除され、別の 2 本が追加されます。対称差の要素数は 4 です。重心は更新前後とも原点であり、エッジ数もともに 2 本です。
\mathbf{V}_0=(0,0),\qquad A_1=\sqrt{4}=2,\qquad \Delta K_1=0
| 確認対象 | 実験 1 の期待値 | 実験 2 の期待値 |
|---|---|---|
| 更新前のエッジ | {(0, 1), (2, 3)} |
{(0, 1), (2, 3)} |
| 更新後のエッジ | {(0, 1), (2, 3)} |
{(0, 2), (1, 3)} |
| 重心速度 | (1.0, 0.0) | (0.0, 0.0) |
| 対称差の要素数 | 0 | 4 |
| 構造振幅 | 0.0 | 2.0 |
| エッジ数の変化 | 0 | 0 |
実験 1 と実験 2 は、観測処理を検証するために指定した座標系列です。各配置は二つの連結成分を持つ固定粒子集合であり、単一の集団が相互作用によって形成された結果を表すものではありません。
相反・非相反相互作用で粒子を更新する
二粒子に働く作用の定義
次に、位置更新を相互作用から計算します。二粒子 A、B に対し、A から B へ向かう単位ベクトルを定義します。
\mathbf{e}_t=\frac{\mathbf{x}_{B,t}-\mathbf{x}_{A,t}}
{\|\mathbf{x}_{B,t}-\mathbf{x}_{A,t}\|_2}
作用範囲内で一定の大きさの力が働く簡略モデルを使用します。粒子 A が B から受ける作用の係数を $g_{AB}$、B が A から受ける作用の係数を $g_{BA}$ とします。正の係数を斥力、負の係数を引力とする規約です[1]。
\mathbf{F}_{A,t}=-g_{AB}\mathbf{e}_t,\qquad
\mathbf{F}_{B,t}=g_{BA}\mathbf{e}_t
作用の和と差は次のようになります。
\begin{aligned}
\mathbf{F}_{A,t}+\mathbf{F}_{B,t}
&=(g_{BA}-g_{AB})\mathbf{e}_t\\
\mathbf{F}_{B,t}-\mathbf{F}_{A,t}
&=(g_{BA}+g_{AB})\mathbf{e}_t
\end{aligned}
各粒子の質量が 1 の場合、合力の半分が重心加速度、作用の差が相対位置の加速度になります。作用範囲内で粒子間距離がゼロでない場合、係数が等しければ合力はゼロになり、係数が異なれば合力は非ゼロになります。
実験 3 では、作用範囲 $r_{\mathrm{int}}$ を 2.0 として step_pair() を経由して pair_forces() に渡します。近接距離 $r_c$ と同じ数値ですが、役割が異なるため別の引数で渡します。粒子間距離が作用範囲以上なら両粒子への作用をゼロとします。距離がゼロの場合は単位ベクトルが未定義となるため、ValueError を発生させます。
速度と位置の離散時間更新
粒子の質量を 1 とし、更新前の位置から計算した作用で速度を更新します。その後、更新した速度を使って位置を更新します。これは半陰的 Euler 法の更新順序です。
\begin{aligned}
\mathbf{v}_{i,t+1}&=\mathbf{v}_{i,t}+\Delta t\mathbf{F}_{i,t}\\
\mathbf{x}_{i,t+1}&=\mathbf{x}_{i,t}+\Delta t\mathbf{v}_{i,t+1}
\end{aligned}
step_pair() は更新前の位置配列と速度配列を直接変更せず、新しい配列を返します。比較する相互作用条件の開始状態を一致させるためです。pair_forces() と step_pair() は二粒子専用で、位置配列と速度配列は形状 (2, 2) の有限値の浮動小数点配列を前提とします。時間刻みと作用範囲には有限の正値を、相互作用係数には有限値を使用します。
陽的 Euler 法のように更新前の速度を位置更新に用いる場合、初期速度がゼロなら最初の位置更新で粒子は移動しません。以下の期待値は、速度を先に更新する上記の式に対応します。
実験 3：相反的な引力と非相反的な運動
初期位置は A が (0, 0)、B が (1, 0) です。両粒子の初期速度をゼロとし、時間刻みを 0.1、作用範囲と近接距離を 2.0 に設定して 1 ステップ更新します。
| 条件 | $g_{AB}$ | $g_{BA}$ |
|---|---|---|
| 相反的な引力 | -1 | -1 |
| 非相反的な運動 | -1 | 1 |
相反的な引力では、A に (1, 0)、B に (-1, 0) の力が働きます。1 ステップ後の期待値は次のとおりです。
\begin{aligned}
\mathbf{v}_{A,1}&=(0.1,0),&\mathbf{v}_{B,1}&=(-0.1,0)\\
\mathbf{x}_{A,1}&=(0.01,0),&\mathbf{x}_{B,1}&=(0.99,0)
\end{aligned}
重心は (0.5, 0) に保たれ、粒子間距離は 0.98 になります。
非相反的な条件では、両粒子に (1, 0) の力が働きます。
\begin{aligned}
\mathbf{v}_{A,1}&=(0.1,0),&\mathbf{v}_{B,1}&=(0.1,0)\\
\mathbf{x}_{A,1}&=(0.01,0),&\mathbf{x}_{B,1}&=(1.01,0)
\end{aligned}
重心は (0.51, 0) へ移動し、粒子間距離は 1.0 に保たれます。
いずれの条件でも粒子間距離は近接距離 2.0 未満です。このため、更新前後のエッジ集合はともに {(0, 1)} となります。
| 確認対象 | 相反的な引力 | 非相反的な運動 |
|---|---|---|
| 重心速度 | (0.0, 0.0) | (0.1, 0.0) |
| 粒子間距離 | 0.98 | 1.00 |
| 構造振幅 | 0.0 | 0.0 |
| エッジ数の変化 | 0 | 0 |
相反的な引力では、粒子間距離が減少しても重心と近接関係は維持されます。非相反的な条件では、近接関係を維持したまま重心が移動します。距離の変化は二値の近接グラフが測定する関係更新とは別の量です。
検証スクリプトと実行結果
スクリプト全体
以下を particle_structure.py として保存します。既知の座標系列と相互作用から計算した座標系列の両方を observe_transition() に渡す構成です。
数値の比較には np.testing.assert_allclose() を使用します。この関数は絶対許容誤差と相対許容誤差による条件で値を比較し、条件を満たさなければ AssertionError を発生させます。strict=True を指定すると形状とデータ型の一致も要求できます[4]。
今回は少数の粒子を使った固定入力の計算なので、rtol=0.0、atol=1e-12 とします。ゼロを期待する値も絶対許容誤差によって比較します。equal_nan=False により NaN を一致として扱いません。集合と整数については完全一致を確認します。
import numpy as np
def proximity_edges(positions, cutoff):
return {
(i, j)
for i in range(len(positions))
for j in range(i + 1, len(positions))
if np.linalg.norm(positions[j] - positions[i]) < cutoff
}
def centroid_velocity(before, after, dt):
if before.shape != after.shape:
raise ValueError("Particle arrays must have the same shape.")
if dt <= 0:
raise ValueError("dt must be positive.")
return (after.mean(axis=0) - before.mean(axis=0)) / dt
def structural_amplitude(before_edges, after_edges):
return np.sqrt(len(before_edges ^ after_edges))
def pair_forces(positions, g_ab, g_ba, interaction_radius):
displacement = positions[1] - positions[0]
distance = np.linalg.norm(displacement)
if distance == 0.0:
raise ValueError("Coincident particles have no unit direction.")
if distance >= interaction_radius:
return np.zeros_like(positions)
direction = displacement / distance
return np.array([-g_ab * direction, g_ba * direction])
def step_pair(positions, velocities, dt, g_ab, g_ba,
interaction_radius):
forces = pair_forces(positions, g_ab, g_ba, interaction_radius)
next_velocities = velocities + dt * forces
next_positions = positions + dt * next_velocities
return next_positions, next_velocities
def observe_transition(before, after, dt, cutoff):
before_edges = proximity_edges(before, cutoff)
after_edges = proximity_edges(after, cutoff)
return {
"velocity": centroid_velocity(before, after, dt),
"before_edges": before_edges,
"after_edges": after_edges,
"changed_edges": len(before_edges ^ after_edges),
"amplitude": structural_amplitude(before_edges, after_edges),
"edge_count_delta": len(after_edges) - len(before_edges),
}
def check_close(actual, expected):
np.testing.assert_allclose(
actual, expected, rtol=0.0, atol=1e-12, equal_nan=False,
strict=True,
)
def check_observation(result, velocity, before_edges, after_edges,
changed_edges, amplitude, edge_count_delta):
check_close(result["velocity"], np.array(velocity, dtype=np.float64))
assert result["before_edges"] == before_edges
assert result["after_edges"] == after_edges
assert result["changed_edges"] == changed_edges
check_close(result["amplitude"], np.float64(amplitude))
assert result["edge_count_delta"] == edge_count_delta
def print_result(name, result):
vx, vy = result["velocity"]
print(
f"{name}: velocity=({vx:.1f}, {vy:.1f}), "
f"amplitude={result['amplitude']:.1f}, "
f"changed_edges={result['changed_edges']}, "
f"edge_count_delta={result['edge_count_delta']}"
)
def initial_positions():
return np.array([
[-0.5, -1.5],
[ 0.5, -1.5],
[-0.5, 1.5],
[ 0.5, 1.5],
], dtype=np.float64)
def test_translation():
before = initial_positions()
after = before + np.array([1.0, 0.0])
result = observe_transition(before, after, 1.0, 2.0)
edges = {(0, 1), (2, 3)}
check_observation(result, (1.0, 0.0), edges, edges, 0, 0.0, 0)
print_result("translation", result)
def test_rearrangement():
before = initial_positions()
after = np.array([
[-1.5, -0.5],
[ 1.5, -0.5],
[-1.5, 0.5],
[ 1.5, 0.5],
], dtype=np.float64)
result = observe_transition(before, after, 1.0, 2.0)
check_observation(
result, (0.0, 0.0), {(0, 1), (2, 3)},
{(0, 2), (1, 3)}, 4, 2.0, 0,
)
print_result("rearrangement", result)
def test_pair_interaction():
positions = np.array([[0.0, 0.0], [1.0, 0.0]], dtype=np.float64)
velocities = np.zeros_like(positions)
initial_x = positions.copy()
initial_v = velocities.copy()
cases = [
("reciprocal", -1.0, [[0.01, 0.0], [0.99, 0.0]],
[[0.1, 0.0], [-0.1, 0.0]], (0.0, 0.0), 0.98),
("nonreciprocal", 1.0, [[0.01, 0.0], [1.01, 0.0]],
[[0.1, 0.0], [0.1, 0.0]], (0.1, 0.0), 1.0),
]
for name, g_ba, expected_x, expected_v, expected_cv, expected_d in cases:
after, next_v = step_pair(positions, velocities, 0.1, -1.0,
g_ba, 2.0)
check_close(after, np.array(expected_x, dtype=np.float64))
check_close(next_v, np.array(expected_v, dtype=np.float64))
np.testing.assert_array_equal(positions, initial_x)
np.testing.assert_array_equal(velocities, initial_v)
distance = np.linalg.norm(after[1] - after[0])
check_close(distance, np.float64(expected_d))
result = observe_transition(positions, after, 0.1, 2.0)
edges = {(0, 1)}
check_observation(result, expected_cv, edges, edges, 0, 0.0, 0)
print_result(name, result)
print(f" distance={distance:.2f}")
def test_boundaries():
at_cutoff = np.array([[0.0, 0.0], [2.0, 0.0]])
outside = np.array([[0.0, 0.0], [3.0, 0.0]])
assert proximity_edges(at_cutoff, 2.0) == set()
for positions in (at_cutoff, outside):
check_close(pair_forces(positions, -1.0, 1.0, 2.0),
np.zeros_like(positions))
with np.testing.assert_raises(ValueError):
pair_forces(np.zeros((2, 2)), -1.0, 1.0, 2.0)
print("boundary checks: PASS")
def main():
test_translation()
test_rearrangement()
test_pair_interaction()
test_boundaries()
print("All checks passed.")
if __name__ == "__main__":
main()
実行方法と期待値の確認
PYTHONOPTIMIZE=0 python particle_structure.py
Python の最適化オプション -O を付けず、環境変数 PYTHONOPTIMIZE=0 を指定して実行します。最適化を無効にし、assert による集合と整数の比較が有効な状態で検証するためです。
実行結果は次のとおりです。
translation: velocity=(1.0, 0.0), amplitude=0.0, changed_edges=0, edge_count_delta=0
rearrangement: velocity=(0.0, 0.0), amplitude=2.0, changed_edges=4, edge_count_delta=0
reciprocal: velocity=(0.0, 0.0), amplitude=0.0, changed_edges=0, edge_count_delta=0
distance=0.98
nonreciprocal: velocity=(0.1, 0.0), amplitude=0.0, changed_edges=0, edge_count_delta=0
distance=1.00
boundary checks: PASS
All checks passed.
三つの実験について、座標、速度、粒子間距離、更新前後のエッジ集合、重心速度、構造振幅、エッジ数の変化を、それぞれ該当する期待値と比較しています。実験 3 では相互作用条件ごとに、入力の位置配列と速度配列が変更されていないことも確認します。
test_boundaries() は、次の境界条件を追加で検証します。
| 確認条件 | 期待値 |
|---|---|
| 粒子間距離が近接距離 2.0 と等しい | エッジ集合は空集合 |
| 粒子間距離が作用範囲 2.0 以上 | 両粒子への作用がゼロ |
| 二粒子が同一座標に存在する |
ValueError が発生 |
実行ではすべての検証を通過しました。表示された構造振幅と重心速度を整理すると次のようになります。
| 実験条件 | 重心速度 | 構造振幅 | エッジ数の変化 |
|---|---|---|---|
| 平行移動 | (1.0, 0.0) | 0.0 | 0 |
| 近接関係の更新 | (0.0, 0.0) | 2.0 | 0 |
| 相反的な引力 | (0.0, 0.0) | 0.0 | 0 |
| 非相反的な運動 | (0.1, 0.0) | 0.0 | 0 |
平行移動と関係更新を比較すると、エッジ数の変化量だけでは近接相手の入れ替わりを検出できないことが確認できます。二粒子の実験では、構造振幅がゼロのままでも、重心速度と粒子間距離を併用することで運動の違いを確認できました。
近接グラフを用いた測定の留意事項
今回の観測処理を別の粒子シミュレーションへ適用するときは、以下の条件を確認します。
- 粒子 ID と配列の行の対応を時刻間で維持します。特定の粒子集合を追跡する場合は、その集合を明示してから測定します。
- 近接距離は生成規則とは独立した観測条件として記録します。閾値付近では浮動小数点誤差によりエッジの判定が変化する場合があるため、比較規則と数値精度を揃えます。
- 重心速度と構造振幅は同じ観測区間に対応付けます。観測間隔が長いと、その区間内で追加された後に削除されたエッジなどの途中経過は記録に残りません。
- 粒子間距離の連続的な変化を評価する場合は、二値グラフに加えて距離そのものを記録します。構造振幅は閾値をまたいだ近接関係の更新を測ります。
- 周期境界を持つモデルでは、境界を考慮した距離計算と重心の計算方法が必要です。今回の座標差と算術平均は通常のユークリッド空間を前提としています。
粒子数が異なる集合間で構造振幅を比較する場合にも注意が必要です。無向エッジ候補の数は $N(N-1)/2$ なので、同じ割合の関係更新であっても構造振幅の絶対値は異なり得ます。候補となる粒子対のうち近接関係が更新された割合を併記する方法があります。
R_{t+1}=\frac{|E_{t+1}\triangle E_t|}{\binom{N}{2}},\qquad N\geq2
この割合はエッジ候補数を分母とする指標です。実際に存在するエッジの本数を分母とする指標とは異なります。
今回の proximity_edges() はすべての粒子対を列挙します。距離計算の回数は $N(N-1)/2$、計算量は $O(N^2)$ です。粒子数を増やす場合は、近接探索の計算量も評価対象にする必要があります。
実験 3 の二条件では、係数の反対称成分だけでなく対称成分も異なります。相互作用の非相反性だけの寄与を検証するには、対称成分を固定して反対称成分を変更する比較が必要です。今回確認した範囲は、指定した初期状態と更新則における 1 ステップの運動と観測値です。持続的な追跡運動や外乱への安定性は、別途時間発展を調べる必要があります。
まとめ
既知の座標系列で検証した観測関数を、相互作用による位置更新にもそのまま適用できました。
粒子シミュレーションを評価する際には、位置更新処理とは独立して、固定粒子 ID に基づく近接グラフの構築と重心速度・構造振幅の計算を実装します。重心速度は平均位置の移動を、構造振幅はエッジの追加・削除を測ります。必要に応じて粒子間距離も記録すれば、二値グラフでは検出できない連続的な配置変化を評価できます。
参考文献
- id774, 集団構造の形成・持続・自己維持は別々の条件で成立する（2026-11-12）. https://blog.id774.net/entry/2026/11/12/5887/
- NumPy Developers, numpy.linalg.norm. https://numpy.org/doc/2.2/reference/generated/numpy.linalg.norm.html
- id774, 構造振動モデルの定義を測定と検証まで含めて整理する（2026-11-11）. https://blog.id774.net/entry/2026/11/11/5869/
- NumPy Developers, numpy.testing.assert_allclose. https://numpy.org/doc/2.2/reference/generated/numpy.testing.assert_allclose.html