diff --git a/.github/workflows/octave_reference.yml b/.github/workflows/octave_reference.yml index 14fddb7..eb4a50b 100644 --- a/.github/workflows/octave_reference.yml +++ b/.github/workflows/octave_reference.yml @@ -41,7 +41,7 @@ jobs: - run: pip install -e '.[test]' # cross-validates PyBMD against the real bmd.m/cbmd.m under Octave; see - # docs/octave_cross_validation.md for the method and measured tables + # tests/octave/octave_cross_validation.md for the method and measured tables - run: pytest tests/test_octave_reference.py -q env: PYBMD_REQUIRE_OCTAVE_REF: '1' diff --git a/.gitignore b/.gitignore index b8c6e74..3bf664d 100644 --- a/.gitignore +++ b/.gitignore @@ -10,10 +10,9 @@ htmlcov/ # results written by fit() and by the examples bmd_results/ example*_out/ -cylinder_sumdiff_out/ # refs/ holds copies of external reference material used only for local/CI -# cross-validation (see docs/octave_cross_validation.md). The MATLAB source +# cross-validation (see tests/octave/octave_cross_validation.md). The MATLAB source # itself (refs/bmd) is a git submodule -- tracked as a pointer, never # vendored -- so its own license (research/non-commercial) never applies to # this MIT-licensed repo. Everything else under refs/ (the paper preprint, diff --git a/CLAUDE.md b/CLAUDE.md index 2a4e279..e4855c5 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -14,13 +14,13 @@ pip install -e '.[mpi,io,test]' # editable install; extras: mpi, io (.mat/. git submodule update --init # populate refs/bmd (the MATLAB reference), only needed # for tests/test_octave_reference.py -- see tests/CLAUDE.md -pytest # full suite, ~90 s, 169 tests: 13 `slow` (Octave - # cross-validation, figure regeneration), 1 `mpi` +pytest # full suite, 143 tests: 12 `slow` (Octave + # cross-validation), 1 `mpi` pytest -m "not slow and not mpi" # fast subset, ~30 s pytest tests/optimizers -q # one directory (numerical-radius solver tests, # one test function per file) -pytest tests/test_bmd_serial.py::test_bispectrum_matches_closed_form -q # one test -pytest -k "conjugate or closed_form" # by name +pytest tests/test_hypothesis.py::test_resonant_triad_is_detected -q # one test +pytest -k "hypothesis or matlab_compat" # by name python -m pyflakes pybmd/ tests/ examples/ # only linter used; ignore the # "f-string is missing placeholders" @@ -38,7 +38,13 @@ BLAS reorders reductions — so set the same when comparing MPI runs by hand. Examples run from any directory (the fixture path is resolved relative to `examples/data.py`): `MPLBACKEND=Agg python examples/example1_cylinder.py`; they write `example*_out/` in the working -directory. +directory. Examples 4 and 5 compare against Schmidt (2020)'s figures: +`examples/example4_hypothesis_testing.py` (surrogate data, Figs. 4-5; its helpers are also +imported by `tests/test_hypothesis.py`) and `examples/example5_cylinder_paper.py` (the +full-resolution cylinder bispectrum and the modes of its labelled triads, Figs. 7-9, ~2-4.5 min, +needs the `refs/bmd` submodule; see `examples/example5_cylinder_paper.md`). The Octave report +figures are regenerated by +`python tests/octave/build_report.py` into `tests/octave/figures/`. ## Architecture diff --git a/README.md b/README.md index f991a40..33c12b6 100644 --- a/README.md +++ b/README.md @@ -9,8 +9,7 @@ associated with it, distinguishing sum- from difference-interactions and produci maps that identify the regions of nonlinear coupling. The architecture follows [PySPOD](https://github.com/MathEXLab/PySPOD): a `params`-dict-driven -`Base`/`Standard` class pair, an optional MPI communicator, disk-backed mode storage, and a YAML -config reader. +`Base`/`Standard` class pair, an optional MPI communicator and disk-backed mode storage. ``` f2 or l @@ -114,7 +113,10 @@ cbmd = Cross(params=dict(params, state_idx=[0], qr_idx=[[1, 2]]), ``` See [`examples/`](examples/) for the three worked cases, which mirror `example1.m`–`example3.m` of -the original MATLAB implementation. +the original MATLAB implementation, and for reproductions of Schmidt (2020)'s figures: +`example4_hypothesis_testing.py` (surrogate data, Figs. 4 and 5) +and `example5_cylinder_paper.py` (cylinder-wake mode bispectrum and modes, Figs. 7-9; see +`example5_cylinder_paper.md`). ## Parameters @@ -183,7 +185,7 @@ is the default. `MengiOverton` bug-for-bug, confirmed live against the real MATLAB source under Octave to a few micro-relative on well-scaled problems. It exists **only** to reproduce a specific published MATLAB result — it reproduces a confirmed -under-estimation bug and should never be used to analyse new data. See `docs/octave_cross_validation.md` +under-estimation bug and should never be used to analyse new data. See `tests/octave/octave_cross_validation.md` for the measured figures and `pybmd.bmd.optimizers.mengi_overton`'s docstring for the caveats. ## Testing @@ -193,12 +195,10 @@ pytest # everything, ~90 s (Octave cross-validation, pytest -m "not slow and not mpi" # fast subset, ~30 s ``` -The suite verifies the bispectrum against a **closed-form analytic result** — for an on-grid, -boxcar-windowed, block-random-phase signal, `L(k,l) = (a_k a_l a_{k+l} / 8) Σ w conj(φ_{k+l}) φ_k φ_l` -for every triad — as well as conjugate symmetry, exact triad counts, CBMD reducing to BMD when the -three variables coincide, bit-identical results between `mpirun -n 1` and `-n 2`, and a -regression against the original MATLAB implementation run live under Octave on the cylinder-wake -dataset (see [`tests/CLAUDE.md`](tests/CLAUDE.md)). +The suite checks the numerical-radius solvers against brute force, reproduces Schmidt (2020)'s +hypothesis test on surrogate data, asserts bit-identical results between `mpirun -n 1` and `-n 2`, +and regresses `L`, `T`, the modes and CBMD against the original MATLAB implementation run live +under Octave on the cylinder-wake dataset (see [`tests/CLAUDE.md`](tests/CLAUDE.md)). ## References diff --git a/docs/bmd_theory_to_implementation.md b/docs/bmd_theory_to_implementation.md new file mode 100644 index 0000000..189a370 --- /dev/null +++ b/docs/bmd_theory_to_implementation.md @@ -0,0 +1,518 @@ +# BMD の理論と PyBMD 実装の対応 + +Schmidt (2020, *Nonlinear Dynamics*) の Bispectral Mode Decomposition (BMD) を、理論の各ステップが +PyBMD のどのコードに対応するかという観点で整理した文書です。式は PyBMD（および MATLAB 版 +`refs/bmd/bmd.m`）が**実際に計算している形**で書いています。論文との表記の違いは §5.3 と §10 にまとめました。 + +実装全体の不変条件や MATLAB 版からの逸脱の詳細は [pybmd/bmd/CLAUDE.md](../pybmd/bmd/CLAUDE.md) を参照してください。 + +--- + +## 0. 全体の流れ + +| 理論のステップ | 数式（概略） | 実装 | +| --- | --- | --- | +| 1. 変動成分をとる | $q' = q - \bar q$ | [base.py:513](../pybmd/bmd/base.py#L513) `select_mean`, [base.py:566](../pybmd/bmd/base.py#L566) | +| 2. ブロック分割（Welch 法） | $N_\mathrm{blk}$ 個の実現 $q^{[i]}$ | [base.py:399-406](../pybmd/bmd/base.py#L399-L406), [base.py:550](../pybmd/bmd/base.py#L550) `_get_block` | +| 3. 窓掛け＋時間 DFT | $\hat q^{[i]}_k$ | [base.py:557](../pybmd/bmd/base.py#L557) `_compute_blocks` | +| 4. トライアド $(k,l,k+l)$ の列挙 | $f_k+f_l=f_{k+l}$ | [utils.py:299](../pybmd/bmd/utils.py#L299) `triad_indices` | +| 5. 実現行列の組み立て | $\hat Q_{k+l},\ \hat Q_{k\circ l}=\hat Q_k\circ\hat Q_l$ | [standard.py:60](../pybmd/bmd/standard.py#L60) `_triad_matrices` | +| 6. バイスペクトル密度行列 | $\mathbf B = \hat Q_{k+l}^H \mathbf W \hat Q_{k\circ l}/N_\mathrm{blk}$ | [base.py:646](../pybmd/bmd/base.py#L646) | +| 7. 数値半径の最大化 | $\mathbf a_1=\arg\max_{\lVert\mathbf a\rVert=1}\lvert\mathbf a^H\mathbf B\mathbf a\rvert$ | [optimizers.py:340](../pybmd/bmd/optimizers.py#L340) `solve` | +| 8. モード・モードバイスペクトル | $\lambda_1,\ \phi_{k+l},\ \phi_{k\circ l}$ | [base.py:656-679](../pybmd/bmd/base.py#L656-L679) | + +呼び出し順は `Base.fit()` → `_initialize` → `_compute_qhat` → `_triad_loop` → `_store_and_save` です +（[standard.py:24](../pybmd/bmd/standard.py#L24)）。 + +### 記号と変数名の対応 + +| 記号 | 意味 | 変数名 | +| --- | --- | --- | +| $N_t$ | スナップショット数 | `nt` | +| $N_\mathrm{fft}$ | 1 ブロックのスナップショット数 | `n_dft` | +| $N_\mathrm{ovlp}$ | ブロックの重なり（スナップショット数） | `n_overlap` | +| $N_\mathrm{blk}$ | ブロック数（＝実現の数） | `n_blocks` | +| $\Delta t$ | 時間刻み | `time_step`, `dt` | +| $n$ | 空間点数 × 変数数（1 モードの長さ） | `nxv = nx * nv` | +| $k, l$ | 符号付き周波数インデックス | `triads.k`, `triads.l`, `triads.kl` | +| $\hat Q_k$ | 周波数 $k$ の実現行列 $(n\times N_\mathrm{blk})$ | `q_hat[row]`, `q1`, `q2`, `q3` | +| $\mathbf W$ | 空間内積の重み（対角） | `weights`（形状 `(n, 1)`） | +| $\mathbf B$ | バイスペクトル密度行列 $(N_\mathrm{blk}\times N_\mathrm{blk})$ | `B` | +| $\mathbf a_1$ | 展開係数（最適化ベクトル） | `a`, `coeffs[i]` | +| $\lambda_1$ | 複素モードバイスペクトル | `r`, `L[f1_idx, f2_idx]` | +| $\phi_{k+l}$ | バイスペクトルモード | `modes[..., 0, ...]` | +| $\phi_{k\circ l}$ | 相互周波数場 (cross-frequency field) | `modes[..., 1, ...]` | + +--- + +## 1. データと前処理 + +### 1.1 データの形 + +データは常に `(nt, *xshape, n_variables)`、つまり **時間が先頭、変数が末尾** です。 +内部では各スナップショットを C order で長さ $n = n_x n_v$ のベクトルに平坦化します +（[base.py:555](../pybmd/bmd/base.py#L555) の `reshape(self._n_dft, -1)`）。 + +### 1.2 平均の除去 + +BMD は変動成分 $q'(\mathbf x,t)=q(\mathbf x,t)-\bar q(\mathbf x)$ の相関を見る手法です。`mean_type` で選びます +（[base.py:513-543](../pybmd/bmd/base.py#L513-L543)）。 + +| `mean_type` | 引く量 | 備考 | +| --- | --- | --- | +| `'longtime'`（既定） | 全時間平均 $\bar q$ | MATLAB 版の既定と同じ | +| `'blockwise'` | 各ブロックの時間平均 | [base.py:568](../pybmd/bmd/base.py#L568) | +| `'zero'`, `'none'` | 何も引かない | 警告が出る | +| 引数 `mean=` | ユーザ指定の平均 | 形状 `(*xshape, nv)` 必須 | + +### 1.3 空間内積の重み $\mathbf W$ + +2 つの場の内積を +$$ +\langle \mathbf u, \mathbf v\rangle_{\mathbf W} = \mathbf u^H \mathbf W \mathbf v = \sum_j w_j\, \overline{u_j}\, v_j +$$ +と定義します。$w_j$ は通常、求積（台形則など）の体積要素です。 + +- 生成: [pybmd/utils/weights.py](../pybmd/utils/weights.py) の `uniform`, `trapz_2d`, `trapz_3d` +- 形状チェック: [base.py:496](../pybmd/bmd/base.py#L496)。平坦ベクトルは受け付けません（並び順の曖昧さでモードが壊れるのを防ぐため） +- 平坦化: [base.py:424](../pybmd/bmd/base.py#L424)（データと同じ C order） + +$\mathbf B$ は $\mathbf W$ に線形なので、重みの規約を変えると $|\lambda_1|$ 全体が定数倍されます +（[examples/example5_cylinder_paper.md](../examples/example5_cylinder_paper.md) の "Spatial weight" 節）。 + +--- + +## 2. スペクトル推定：ブロック分割・窓・DFT + +### 2.1 ブロック分割 + +期待値 $E[\cdot]$ を、時系列を重なりのあるブロックに分けたアンサンブル平均で近似します（Welch 法）。 + +$$ +N_\mathrm{blk} = \left\lfloor \frac{N_t - N_\mathrm{ovlp}}{N_\mathrm{fft} - N_\mathrm{ovlp}} \right\rfloor +$$ + +- 実装: [base.py:400-402](../pybmd/bmd/base.py#L400-L402)。$N_\mathrm{blk}<2$ ならエラーになります。 +- `overlap` は**パーセント**（既定 50）、`n_overlap` は**スナップショット数**で、後者が優先されます + （[base.py:359](../pybmd/bmd/base.py#L359)）。$N_\mathrm{ovlp}=\lfloor N_\mathrm{fft}\cdot\text{overlap}/100\rfloor$ です。 +- $i$ 番目のブロックの開始位置は $\min(iN_\mathrm{step}+N_\mathrm{fft},N_t)-N_\mathrm{fft}$ です（$N_\mathrm{step}=N_\mathrm{fft}-N_\mathrm{ovlp}$）。 + 最後のブロックがデータの末尾に揃えられます（[base.py:552](../pybmd/bmd/base.py#L552)）。 + +### 2.2 窓と DFT + +ブロック $i$ の $j$ 番目のスナップショットに窓 $w_j$ を掛けて DFT をとります。 + +$$ +\hat q^{[i]}_k = \frac{1}{\bar w\,N_\mathrm{fft}} \sum_{j=0}^{N_\mathrm{fft}-1} w_j\, q'^{[i]}_j\, e^{-2\pi \mathrm i\, jk/N_\mathrm{fft}}, +\qquad \bar w = \frac{1}{N_\mathrm{fft}}\sum_j w_j +$$ + +- $1/\bar w$ は窓による振幅低下の補正です（`win_weight`、[base.py:428](../pybmd/bmd/base.py#L428)）。 + この規格化により、周波数 $f_k$ で振幅 $A$ の正弦波（窓がその周波数に合っている場合）は $|\hat q_k| = A/2$ になります。 +- DFT と `fftshift`: [base.py:581-582](../pybmd/bmd/base.py#L581-L582) +- 窓: `'hamming'`（既定）/`'hann'`/`'boxcar'`/任意の配列（[utils.py:76](../pybmd/bmd/utils.py#L76)） + +**両側スペクトルが必須です。** 差の相互作用（$l<0$）は負の周波数を使うので、`rfft` は使いません。 + +### 2.3 周波数軸 + +`fftshift` 後の行番号と符号付きインデックス $k$、物理周波数の関係は次のとおりです（[utils.py:105](../pybmd/bmd/utils.py#L105) `freq_axis`）。 + +$$ +k \in \{-N_\mathrm{fft}/2, \dots, N_\mathrm{fft}/2-1\}\ (\text{偶数 } N_\mathrm{fft}),\qquad +f_k = \frac{k}{N_\mathrm{fft}\Delta t},\qquad k_\mathrm{Nyq} = N_\mathrm{fft}/2 +$$ + +| 量 | 変数 | +| --- | --- | +| 符号付きインデックス $k$（行ごと） | `triads.f_idx` | +| 物理周波数 $f_k$ | `triads.freq`（`bmd.freq`） | +| Nyquist インデックス | `triads.f_nyq_idx` | +| $k$ → 行番号 | `triads.row_of(k)` | + +`dt` を省略すると $\Delta t = 1/N_\mathrm{fft}$ とされ、$f_k = k$ になります（MATLAB 版と同じ）。 + +### 2.4 実現行列 $\hat Q_k$ + +周波数 $k$ の全ブロックの DFT を列に並べたものです。 + +$$ +\hat Q_k = \begin{bmatrix} \hat q^{[1]}_k & \hat q^{[2]}_k & \cdots & \hat q^{[N_\mathrm{blk}]}_k \end{bmatrix} \in \mathbb C^{n\times N_\mathrm{blk}} +$$ + +実装では dict `q_hat[row]` に格納されます（[base.py:584](../pybmd/bmd/base.py#L584) `_compute_qhat`）。 +メモリ節約のため、どれかのトライアドが参照する行（`triads.freq_needed`）だけを保持します。 + +--- + +## 3. トライアドと $f_1$–$f_2$ 平面 + +### 3.1 トライアド条件 + +2 次の非線形項（例：$\mathbf u\cdot\nabla\mathbf u$）は、周波数 $f_k$ と $f_l$ の成分の積から周波数 $f_k+f_l$ の成分を作ります。 +このため BMD は、周波数の三つ組 + +$$ +(f_k,\ f_l,\ f_{k+l}),\qquad f_k + f_l - f_{k+l} = 0 +$$ + +（**トライアド**）を単位として解析します。インデックスで書けば $(k, l, k+l)$ です。 +$l<0$ なら $f_{k+l} = f_k - |f_l|$ なので差の相互作用を表します。 + +### 3.2 計算するトライアドの選択 + +[utils.py:299](../pybmd/bmd/utils.py#L299) `triad_indices` は、次の条件をすべて満たす $(k,l)$ を列挙します（[utils.py:331-333](../pybmd/bmd/utils.py#L331-L333)）。 + +- $|k+l| < k_\mathrm{Nyq}$（和の周波数も軸上にあること） +- $|k|, |l| \le$ `max_freq_idx`（既定は Nyquist まで） +- $(k,l)$ が指定した `regions` のいずれかに入っていること + +``` + f2 or l + ^ + ________| + |\ |\ + | \ 7 | \ + | 6 \ | 8 /\ + | \| / 1 \ + ----+-------+-------+-> f1 or k + \ 5 / |\ | + \/ 4 | \ 2 | + \ | 3 \ | + \|______\| +``` + +各領域の判定式は [utils.py:126](../pybmd/bmd/utils.py#L126) `_region_masks` にあり、MATLAB 版と同一です。 +`regions` は **1 始まり**（図のラベル番号）です。 + +### 3.3 既定値が `regions=[1, 2]` である理由（対称性） + +$\mathbf B$ は次の 2 つの対称性を持ちます。 + +1. **$k \leftrightarrow l$ の入れ替え**：$\hat Q_k\circ\hat Q_l = \hat Q_l\circ\hat Q_k$ なので $\mathbf B(k,l)=\mathbf B(l,k)$。したがって $\lambda_1(k,l)=\lambda_1(l,k)$ です。 +2. **実数データの共役対称性**：$\hat q_{-k} = \overline{\hat q_k}$ より $\mathbf B(-k,-l) = \overline{\mathbf B(k,l)}$。したがって $\lambda_1(-k,-l)=\overline{\lambda_1(k,l)}$ です。 + +この 2 つを使うと、実数データでは領域 1（和の相互作用 $k\ge l\ge 0$）と領域 2（差の相互作用 $k\ge|l|,\ l\le 0$）で平面全体を代表できます。 +**複素数データ**では 2. が成り立たないので、必要に応じて他の領域も指定してください。 + +（乱数データで $L(3,2)=L(2,3)$ と $L(-3,-2)=\overline{L(3,2)}$ を数値的に確認済みです。） + +### 3.4 `Triads` オブジェクト + +| フィールド | 意味 | +| --- | --- | +| `f1_idx`, `f2_idx`, `f3_idx` | トライアド $i$ の $k,\ l,\ k+l$ の**行番号** | +| `k`, `l`, `kl` | 同じく**符号付きインデックス** | +| `f1`, `f2`, `f3` | 同じく物理周波数 | +| `region` | 所属領域（境界で重なる場合は番号の大きい方） | +| `triad_map` | `(n_freq, n_freq)` のトライアド番号表（対象外は −1） | +| `find(k, l)` | $(k,l)$ → トライアド番号 | + +--- + +## 4. 古典バイスペクトルから BMD へ + +### 4.1 古典バイスペクトル + +1 点の信号 $q(t)$ のバイスペクトルは 3 次の相関 + +$$ +S(f_k, f_l) = E\!\left[\hat q_k\, \hat q_l\, \overline{\hat q_{k+l}}\right] +$$ + +です。$f_k, f_l, f_{k+l}$ の位相が $\theta_k+\theta_l-\theta_{k+l}=\text{const}$ で結合している（**2 次の位相結合**）ときだけ、 +ブロック平均で打ち消されずに大きな値が残ります。位相がランダムな成分は平均で消えます。 + +### 4.2 2 次項 $\hat q_{k\circ l}$ + +BMD では、周波数 $k$ と $l$ の成分の各点ごとの積（アダマール積） + +$$ +\hat q_{k\circ l}(\mathbf x) = \hat q_k(\mathbf x)\,\hat q_l(\mathbf x),\qquad \hat Q_{k\circ l} = \hat Q_k \circ \hat Q_l +$$ + +を 2 次の非線形項の代わりに使います。$\hat q_{k+l}$ と $\hat q_{k\circ l}$ の相関は、空間分布まで含めた古典バイスペクトルの一般化になっています。 + +- 実装: [standard.py:73-76](../pybmd/bmd/standard.py#L73-L76)（`q1 * q2` が $\hat Q_{k\circ l}$、`q3` が $\hat Q_{k+l}$） + +### 4.3 バイスペクトル密度行列 $\mathbf B$ + +$$ +\boxed{\ \mathbf B = \frac{1}{N_\mathrm{blk}}\, \hat Q_{k+l}^H\, \mathbf W\, \hat Q_{k\circ l}\ } \in \mathbb C^{N_\mathrm{blk}\times N_\mathrm{blk}}, +\qquad +B_{ij} = \frac{1}{N_\mathrm{blk}} \left\langle \hat q^{[i]}_{k+l},\ \hat q^{[j]}_{k\circ l} \right\rangle_{\mathbf W} +$$ + +- 実装: [base.py:646](../pybmd/bmd/base.py#L646) `B = q_sum.conj().T @ (q_prod * weights) / self._n_blocks` +- MATLAB 版の `B = Q_hat_f3'*bsxfun(@times,Q_hat_f1.*Q_hat_f2,weight)/nBlks` と同じ式です（`'` は共役転置）。 + +**古典バイスペクトルとの関係.** $\mathbf B$ の対角成分は同じブロック内の相関、非対角成分は異なるブロック間の相関です。 + +$$ +\operatorname{tr}\mathbf B = \frac{1}{N_\mathrm{blk}}\sum_i \left\langle \hat q^{[i]}_{k+l},\ \hat q^{[i]}_{k}\circ\hat q^{[i]}_{l}\right\rangle_{\mathbf W} += \sum_j w_j\ \widehat{E}\!\left[\overline{\hat q_{k+l}}\,\hat q_k\hat q_l\right](\mathbf x_j) +$$ + +つまり $\operatorname{tr}\mathbf B$ は、各点の古典バイスペクトル $S(f_k,f_l)$ の推定値を空間積分したものです。 +これは $\mathbf B$ の対角成分（同じブロック内の相関）だけを等しく足した量です。 +BMD は非対角成分（異なるブロック間の相関）も含めて、ブロックの線形結合の係数 $\mathbf a$ を最適化します（次節）。 + +--- + +## 5. 最適化問題：数値半径 + +### 5.1 定式化 + +同じ係数 $\mathbf a\in\mathbb C^{N_\mathrm{blk}}$ で 2 つの場を実現の線形結合として作ります。 + +$$ +\boldsymbol\psi_{k+l} = \hat Q_{k+l}\,\mathbf a,\qquad \boldsymbol\psi_{k\circ l} = \hat Q_{k\circ l}\,\mathbf a +$$ + +このとき + +$$ +\mathbf a^H \mathbf B\, \mathbf a = \frac{1}{N_\mathrm{blk}}\left\langle \boldsymbol\psi_{k+l},\ \boldsymbol\psi_{k\circ l}\right\rangle_{\mathbf W} +$$ + +です。この相関の大きさを最大にする $\mathbf a$ を求めます。 + +$$ +\mathbf a_1 = \arg\max_{\lVert\mathbf a\rVert_2=1} \left|\mathbf a^H \mathbf B\,\mathbf a\right|, +\qquad +\lambda_1 = \mathbf a_1^H \mathbf B\,\mathbf a_1, +\qquad +|\lambda_1| = r(\mathbf B) +$$ + +$r(\mathbf B)$ は $\mathbf B$ の**数値半径**（field of values の最大絶対値）です。$\mathbf B$ はエルミートではないので、固有値問題ではなく数値半径の最大化問題になります。 +制約 $\lVert\mathbf a\rVert_2=1$ は $\mathbf W$ を含まない通常のユークリッドノルムです。 + +- 実装: [base.py:652-655](../pybmd/bmd/base.py#L652-L655) `r, a = optimizers.solve(B, ...)` +- `r` は**複素数** $\lambda_1$ です（絶対値ではありません）。`L` にはこの複素数がそのまま入り、図にするとき $|L|$ をとります（[postproc.py:271](../pybmd/bmd/postproc.py#L271)）。 + +### 5.2 数値半径の性質とソルバ + +回転させたエルミート部分 + +$$ +\mathbf H(\theta) = \tfrac12\left(e^{\mathrm i\theta}\mathbf B + e^{-\mathrm i\theta}\mathbf B^H\right) +$$ + +を使うと、 + +$$ +r(\mathbf B) = \max_{\theta\in[0,2\pi)} \lambda_{\max}\big(\mathbf H(\theta)\big) +$$ + +です。最大を与える $\theta^\ast$ での $\mathbf H(\theta^\ast)$ の最大固有ベクトルが $\mathbf a_1$ になります。 + +| 関数 | 理論上の役割 | +| --- | --- | +| [optimizers.py:51](../pybmd/bmd/optimizers.py#L51) `max_fov(A, theta)` | $\lambda_{\max}(\mathbf H(\theta))$ | +| [optimizers.py:95](../pybmd/bmd/optimizers.py#L95) `_dominant_eigvec(A, phi)` | $\mathbf H(\phi)$ の最大固有ベクトル $\mathbf a$ と $\mathbf a^H\mathbf B\mathbf a$ | +| [optimizers.py:229](../pybmd/bmd/optimizers.py#L229) `mengi_overton` | Mengi & Overton (2005) のレベルセット法。**既定**、大域収束 | +| [optimizers.py:177](../pybmd/bmd/optimizers.py#L177) `simple_iteration` | Watson の単純反復（論文付録 Algorithm 1）。局所解のみ | +| [optimizers.py:142](../pybmd/bmd/optimizers.py#L142) `default_start` | 初期ベクトル（$\theta$ の粗い走査。乱数を使わない） | +| [optimizers.py:107](../pybmd/bmd/optimizers.py#L107) `_pow2_scale` | $\lVert\mathbf B\rVert_1\in(1/2,1]$ への 2 のべき乗スケーリング | + +**Mengi–Overton 法の要点.** レベル $w$ に対し、$\lambda_{\max}(\mathbf H(\theta)) = w$ となる角度 $\theta$ は、一般化固有値問題 + +$$ +\mathbf R(w)\,\mathbf v = \mu\,\mathbf S\,\mathbf v,\qquad +\mathbf R(w)=\begin{bmatrix}2w\mathbf I & -\mathbf B^H\\ \mathbf I & \mathbf 0\end{bmatrix},\quad +\mathbf S=\begin{bmatrix}\mathbf B & \mathbf 0\\ \mathbf 0 & \mathbf I\end{bmatrix} +$$ + +の単位円上の固有値 $\mu = e^{\mathrm i\theta}$ として得られます（[optimizers.py:303-307](../pybmd/bmd/optimizers.py#L303-L307)）。 +交差角で区切られた区間の中点のうち、$w$ を超えるものを次の候補にします。候補がなくなれば、現在のレベルが大域最大です。 + +**Watson の単純反復.** 次の更新を収束するまで繰り返します（[optimizers.py:212-214](../pybmd/bmd/optimizers.py#L212-L214)）。 + +$$ +w_{m} = \mathbf a_m^H\mathbf B\mathbf a_m,\qquad +\mathbf a_{m+1} \propto w_m\,\mathbf B^H\mathbf a_m + \overline{w_m}\,\mathbf B\,\mathbf a_m +$$ + +**スケーリング.** 数値半径は $r(c\mathbf B)=c\,r(\mathbf B)$（$c>0$）を満たし、最大化ベクトルは変わりません。 +実際の $\mathbf B$ は $1/N_\mathrm{blk}$ と重みのため非常に小さく（$\lVert\mathbf B\rVert_1\sim10^{-6}$ など）、MATLAB 版の絶対許容誤差では交差角が全て棄却されて過小評価が起きます。 +PyBMD は 2 のべき乗で正規化してから解き（2 進浮動小数点で誤差なし）、最後に元の $\mathbf B$ で $\mathbf a^H\mathbf B\mathbf a$ を評価し直します。 +この修正を含む MATLAB 版からの逸脱の詳細は [pybmd/bmd/CLAUDE.md](../pybmd/bmd/CLAUDE.md) の "Deviations" 節にあります。 + +`solver='MengiOvertonMATLAB'` は、MATLAB 版の結果（過小評価を含む）を再現したいときだけ使う互換モードです。 + +### 5.3 論文表記との関係（$\mathbf B$ と $\mathbf B^H$） + +論文本文は本リポジトリに含まれていないため、論文中の $\mathbf B$ の積の順序とは照合していません。 +仮に論文が逆順 $\tilde{\mathbf B} = \hat Q_{k\circ l}^H\mathbf W\hat Q_{k+l}/N_\mathrm{blk} = \mathbf B^H$ で定義していても、 + +$$ +\mathbf a^H\mathbf B^H\mathbf a = \overline{\mathbf a^H\mathbf B\,\mathbf a} +$$ + +なので、$|\lambda_1|$（モードバイスペクトル）と $\mathbf a_1$ は同じで、**複素数 $\lambda_1$ の位相だけが共役**になります。 +PyBMD の `L` は MATLAB 版 `bmd.m` と同じ規約（上の $\mathbf B$）です。 + +### 5.4 $\mathbf a$ の位相の不定性 + +$\mathbf a\to e^{\mathrm i\alpha}\mathbf a$ としても $\mathbf a^H\mathbf B\mathbf a$ は変わりません。したがって $\lambda_1$ は一意ですが、 +$\mathbf a_1$ とモードは**単位複素数倍の不定性**を持ちます。モード同士を比べるときは要素ごとではなく +$|\langle\mathbf u,\mathbf v\rangle|/(\lVert\mathbf u\rVert\lVert\mathbf v\rVert)\approx1$ で比べてください。 +また、$\lambda_1$ が縮退していると（複数の $\mathbf a$ が同じ値を与えると）モードは一意に定まりません。 + +--- + +## 6. 出力される量 + +### 6.1 モードバイスペクトル $\lambda_1(f_k, f_l)$ + +トライアド $i$ の $\lambda_1$ は `L[f1_idx[i], f2_idx[i]]` に入ります（[base.py:658](../pybmd/bmd/base.py#L658)）。 +`L` は `(n_freq, n_freq)` の複素配列で、計算しなかった $(k,l)$ は NaN です（[base.py:697-699](../pybmd/bmd/base.py#L697-L699)）。 + +- $|\lambda_1|$ が大きい ⇔ $(k,l,k+l)$ の間に空間的にコヒーレントな 2 次の位相結合がある +- 図: `plot_mode_bispectrum`（既定で $\log|\lambda_1|$） + +### 6.2 モード + +$$ +\phi_{k+l} = \frac{\hat Q_{k+l}\mathbf a_1}{\lVert\hat Q_{k+l}\mathbf a_1\rVert_{\mathbf W}},\qquad +\phi_{k\circ l} = \frac{\hat Q_{k\circ l}\mathbf a_1}{\lVert\hat Q_{k\circ l}\mathbf a_1\rVert_{\mathbf W}},\qquad +\lVert\mathbf u\rVert_{\mathbf W}=\sqrt{\mathbf u^H\mathbf W\mathbf u} +$$ + +| モード | 名称 | インデックス | 実装 | +| --- | --- | --- | --- | +| $\phi_{k+l}$ | バイスペクトルモード（和の周波数の構造） | 0 | [base.py:656](../pybmd/bmd/base.py#L656), [base.py:667](../pybmd/bmd/base.py#L667) | +| $\phi_{k\circ l}$ | 相互周波数場（$k$ と $l$ の積が作る構造） | 1 | [base.py:657](../pybmd/bmd/base.py#L657), [base.py:668](../pybmd/bmd/base.py#L668) | +| $\phi_k$ | 構成モード（PyBMD 独自、`constituent_modes=True`） | 2 | [base.py:669-674](../pybmd/bmd/base.py#L669-L674) | +| $\phi_l$ | 同上 | 3 | 同上 | + +- 正規化: [utils.py:363](../pybmd/bmd/utils.py#L363) `normalize_mode` +- 平坦ベクトル → 場の形 `(*xshape, nv)`: [base.py:190](../pybmd/bmd/base.py#L190) `_unflatten_modes` +- 取得: `bmd.get_modes_at_freqs(k, l)`、`bmd.get_modes_at_triad(i)` +- 図: `plot_triad_modes` は $\mathrm{Re}\,\phi$ と、相互作用の局在を示す $|\phi_{k\circ l}\cdot\phi_{k+l}|$（各点の積の絶対値）を描きます（[postproc.py:417](../pybmd/bmd/postproc.py#L417)）。 + +**注意**: $\phi_{k\circ l} \ne \phi_k\circ\phi_l$ です。$(\hat Q_k\circ\hat Q_l)\mathbf a \ne (\hat Q_k\mathbf a)\circ(\hat Q_l\mathbf a)$ だからです。 +$\phi_k,\phi_l$ は同じ $\mathbf a_1$ で作った別の情報で、$\phi_{k\circ l}$ の分解ではありません。 + +### 6.3 エネルギー輸送項 $T$ + +$$ +T(f_k,f_l) = \frac{1}{N_\mathrm{blk}}\,\mathrm{Re}\left[\left(\hat Q_{k+l}\mathbf a_1\right)^H\left(\hat Q_{k\circ l}\mathbf a_1\right)\right] +$$ + +- 実装: [base.py:662-663](../pybmd/bmd/base.py#L662-L663)。正規化前の $\boldsymbol\psi$ を使います。 +- **重み $\mathbf W$ を含みません。** MATLAB 版と同じで、意図的です。 +- したがって一様重み（$\mathbf W=\mathbf I$）なら $T = \mathrm{Re}\,\lambda_1$ です（数値的に確認済み）。 + 重みがある場合は $\mathbf B$ から $\mathbf W$ を除いた Rayleigh 商の実部になります。 +- 符号付きの量で、$f_{k+l}$ への（正）/からの（負）正味のエネルギー輸送を表します。図は `plot_energy_transfer`。 + +### 6.4 展開係数 $\mathbf a_1$ + +`coeffs` は `(n_triads, n_blocks)` で、各トライアドの $\mathbf a_1$ です（[base.py:664](../pybmd/bmd/base.py#L664)）。 +モードは $\hat Q\,\mathbf a_1$ で再構成できるので、大規模な計算では `save_modes=False` として `coeffs.npy` だけを残せます +（ただし再構成用のヘルパーはまだありません）。 + +### 6.5 保存ファイル + +`savedir/nfft{N}_novlp{M}_nblks{B}/` に保存されます（[base.py:752](../pybmd/bmd/base.py#L752) `_store_and_save`）。 + +| ファイル | 中身 | +| --- | --- | +| `bispectrum.npz` | `L`（$\lambda_1$）, `T`, `freq`, `f_idx` | +| `triads.npz` | `Triads` の全フィールド | +| `coeffs.npy` | $\mathbf a_1$ | +| `weights.npy` | 平坦化された $\mathbf W$ | +| `ltm_modes.npy` | 時間平均 $\bar q$ | +| `modes/triad_idx_XXXXXXXX.npy` | トライアドごとのモード `(n_comp, *xshape, nv)` | +| `params_modes.yaml` | パラメータ | + +保存結果の読み込みは [pybmd/bmd/postproc.py](../pybmd/bmd/postproc.py) の `load_results` です。 + +--- + +## 7. Cross-BMD（CBMD） + +### 7.1 理論 + +実際の非線形項は、異なる変数の積の和です（例：$u\,\partial_x u + v\,\partial_y u$）。CBMD は、状態 $s$ と、それを作る積 $q\,r$ の和との相関を見ます。 + +$$ +\hat Q_{k+l} \to \hat S_{k+l},\qquad +\hat Q_{k\circ l} \to \sum_{m} \hat Q^{(m)}_k \circ \hat R^{(m)}_l +$$ + +複数の状態 $s^{(1)},\dots,s^{(n_s)}$ を扱う場合は、それらを縦に積んで 1 本のベクトルにします。 + +$$ +\hat Q_\mathrm{s} = \begin{bmatrix}\hat S^{(1)}_{k+l}\\ \vdots\\ \hat S^{(n_s)}_{k+l}\end{bmatrix},\qquad +\hat Q_\mathrm{qr} = \begin{bmatrix}\sum_m \hat Q^{(1,m)}_k\circ\hat R^{(1,m)}_l\\ \vdots\end{bmatrix},\qquad +\mathbf B = \frac{1}{N_\mathrm{blk}}\hat Q_\mathrm s^H\mathbf W\hat Q_\mathrm{qr} +$$ + +以降（数値半径、モード、$T$）は BMD と同じです。 + +### 7.2 実装 + +| 理論 | 実装 | +| --- | --- | +| $s$ のインデックス | `params['state_idx']`（**0 始まり**、既定 `[0]`） | +| $(q,r)$ の組 | `params['qr_idx']`（列が $q,r,q,r,\dots$ と交互、状態ごとに 1 行。既定 `[[1, 2]]`） | +| $\hat Q_\mathrm s$, $\hat Q_\mathrm{qr}$ の組み立て | [cross.py:153](../pybmd/bmd/cross.py#L153) `Cross._triad_matrices`（和は [cross.py:177-179](../pybmd/bmd/cross.py#L177-L179)） | +| 重み（空間のみ、状態数だけタイル） | [cross.py:128](../pybmd/bmd/cross.py#L128) | +| モードの形 `(*xshape, n_state)` | [cross.py:99](../pybmd/bmd/cross.py#L99) `_unflatten_modes` | + +平坦軸は**状態が最も遅い添字**（`flat = j*nx + p`）です。そのため `_unflatten_modes` は `(n_state, *xshape)` に戻してから状態軸を末尾へ移します。 +CBMD では `normalize_weights` と `constituent_modes` は使えません。 + +--- + +## 8. PyBMD 独自の拡張（MATLAB 版にないもの） + +| 機能 | 理論上の意味 | 実装 | +| --- | --- | --- | +| `constituent_modes=True` | $\phi_k=\hat Q_k\mathbf a_1$, $\phi_l=\hat Q_l\mathbf a_1$ も出力 | [base.py:669-674](../pybmd/bmd/base.py#L669-L674) | +| `normalize_weights=True` | 変数ごとに $w\leftarrow w/\operatorname{var}(q_v)$（異なる単位の変数を揃える） | [weights.py:94](../pybmd/utils/weights.py#L94) | +| `normalize_data=True` | 各ブロック・各点・各変数を標準偏差で割る | [base.py:571-577](../pybmd/bmd/base.py#L571-L577) | +| `mean_type='blockwise'` | ブロックごとの平均を除去 | [base.py:568](../pybmd/bmd/base.py#L568) | +| `window='hann'/'boxcar'` | 窓の選択 | [utils.py:76](../pybmd/bmd/utils.py#L76) | +| MPI 並列 | トライアドをラウンドロビンで分配し `allreduce` | [base.py:638](../pybmd/bmd/base.py#L638), [base.py:690-694](../pybmd/bmd/base.py#L690-L694) | + +`normalize_data=True` は複素数の入力データに対して分散の計算が誤っています（$|x|^2$ ではなく $x^2$ を使っている）。 +実数データには影響しません。詳細は [docs/complex-conjugation-audit.md](complex-conjugation-audit.md)。 + +--- + +## 9. 最小の使用例と理論量の対応 + +```python +from pybmd.bmd.standard import Standard +import pybmd.utils.weights as W + +params = dict(n_dft=64, # N_fft + time_step=dt, # Δt + n_space_dims=2, n_variables=1, + overlap=50, # N_ovlp = 32 + regions=[1, 2], # 和と差の相互作用 + max_freq_idx=12) # |k|,|l| ≤ 12 +bmd = Standard(params, weights=W.trapz_2d(x, y, n_vars=1)).fit(data) # data: (nt, nx, ny, 1) + +i = bmd.find_triad(12, 12) # トライアド (12, 12, 24) +lam = bmd.L[bmd.triads.f1_idx[i], bmd.triads.f2_idx[i]] # λ_1（複素数） +phi = bmd.get_modes_at_triad(i) # phi[0] = φ_{k+l}, phi[1] = φ_{k∘l} +a1 = bmd.coeffs[i] # a_1 +``` + +実例は [examples/example1_cylinder.py](../examples/example1_cylinder.py)（BMD）と [examples/example3_cbmd.py](../examples/example3_cbmd.py)（CBMD）を参照してください。 + +--- + +## 10. 読むときの注意点のまとめ + +- `L` は**複素数** $\lambda_1$。モードバイスペクトルはその絶対値 $|\lambda_1|$。位相の規約は MATLAB 版と同じ（§5.3）。 +- `T` は重みを含まない。一様重みなら $T=\mathrm{Re}\,\lambda_1$。 +- モードは単位複素数倍の不定性を持つ（§5.4）。 +- 実数データなら `regions=[1,2]` で平面全体を代表できる。複素数データでは不十分（§3.3）。 +- `regions` は 1 始まり、`state_idx`/`qr_idx` は 0 始まり。 +- $|\lambda_1|$ の絶対値は重み $\mathbf W$ の規約に比例して変わる。MATLAB 版の図と比べるときは一様重みを使う。 +- 既定ソルバ `MengiOverton` は MATLAB 版の過小評価を修正している。MATLAB 版との数値の差はこれが主因（[tests/octave/octave_cross_validation.md](../tests/octave/octave_cross_validation.md)）。 diff --git a/docs/build_cylinder_figure.py b/docs/build_cylinder_figure.py deleted file mode 100644 index 9a98b7e..0000000 --- a/docs/build_cylinder_figure.py +++ /dev/null @@ -1,234 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -''' -Reproduce ``refs/figures/cylinder_bispectrum_sumdiff.pdf`` (Schmidt 2020, -*Nonlinear Dynamics*) with PyBMD on the cylinder-wake dataset. - -Requires the ``refs/bmd`` submodule populated (``git submodule update --init``) -for ``refs/bmd/wake_Re500.mat``; run from anywhere, paths are resolved -relative to this file. - - MPLBACKEND=Agg python docs/build_cylinder_figure.py - MPLBACKEND=Agg mpirun -n 4 python docs/build_cylinder_figure.py --mpi - -See docs/cylinder_bispectrum.md for how the parameters below were derived -from the published figure. -''' -import argparse -import os -import sys - -import matplotlib -matplotlib.use('Agg') -import matplotlib.pyplot as plt -from matplotlib.cm import ScalarMappable -from matplotlib.colors import LinearSegmentedColormap, Normalize -import numpy as np - -DOCS_DIR = os.path.dirname(os.path.realpath(__file__)) -REPO_ROOT = os.path.realpath(os.path.join(DOCS_DIR, '..')) -FIG_DIR = os.path.join(DOCS_DIR, 'figures', 'cylinder') -sys.path.insert(0, REPO_ROOT) - -from pybmd.bmd.standard import Standard -from pybmd.bmd.postproc import load_results, top_triads -import pybmd.utils.weights as utils_weights -from pybmd.utils.io import read_data - -DEFAULT_DATA = os.path.join(REPO_ROOT, 'refs', 'bmd', 'wake_Re500.mat') - -# the six triads the reference figure circles and labels in panel (b) -TRIADS = ((12, 12), (12, 0), (24, 12), (24, 24), (36, 12), (36, 24)) -# per-triad label placement: (dx, dy) offset in points and horizontal -# alignment. The two triads sharing a row (k=24 and k=36) point their labels -# away from each other so the text doesn't collide in the gap between them; -# (12,12) and (12,0) are pushed right, clear of the y-axis and the f2=0 line -LABEL_LAYOUT = { - (12, 12): ((6, 2), 'left'), (12, 0): ((8, -2), 'left'), - (24, 12): ((-4, 8), 'right'), (24, 24): ((-4, 8), 'right'), - (36, 12): ((4, 8), 'left'), (36, 24): ((4, 8), 'left'), -} - -# sampled from refs/figures/cylinder_bispectrum_sumdiff.pdf's own colorbar, -# rasterized at 400 dpi: jet, with its low (dark blue) end replaced by a fade -# to white so the featureless background reads as blank rather than "cold" -_CBAR_STOPS = ( - (0.000, (1.000, 1.000, 1.000)), (0.074, (0.835, 0.835, 1.000)), - (0.152, (0.482, 0.482, 1.000)), (0.230, (0.075, 0.075, 1.000)), - (0.307, (0.000, 0.345, 1.000)), (0.381, (0.004, 0.733, 1.000)), - (0.459, (0.051, 1.000, 0.953)), (0.537, (0.306, 1.000, 0.694)), - (0.615, (0.694, 1.000, 0.306)), (0.689, (1.000, 1.000, 0.000)), - (0.767, (1.000, 0.545, 0.000)), (0.844, (1.000, 0.090, 0.000)), - (0.922, (0.776, 0.000, 0.000)), (1.000, (0.502, 0.000, 0.000)), -) - - -def _check_prereqs(data_path): - if not os.path.exists(data_path): - sys.exit( - f'Cannot build the figure:\n {data_path} not found -- run ' - '`git submodule update --init` to populate refs/bmd.') - - -def _cylinder_colormap(): - return LinearSegmentedColormap.from_list('cylinder_jet', _CBAR_STOPS) - - -def _parse_args(): - p = argparse.ArgumentParser(description=__doc__) - p.add_argument('--n-dft', type=int, default=480) - p.add_argument('--overlap', type=float, default=50, - help='percent overlap between DFT blocks (default: 50, ' - 'matching the reference)') - p.add_argument('--data', default=DEFAULT_DATA) - p.add_argument('--savedir', default=os.path.join( - REPO_ROOT, 'cylinder_sumdiff_out')) - p.add_argument('--out', default=os.path.join( - FIG_DIR, 'cylinder_bispectrum_sumdiff.png')) - p.add_argument('--dpi', type=int, default=200) - p.add_argument('--clim', type=float, nargs=2, default=None, - metavar=('VMIN', 'VMAX'), - help='override the shared log|lambda_1| colour limits') - p.add_argument('--reuse', action='store_true', - help='load an existing --savedir instead of re-fitting') - p.add_argument('--mpi', action='store_true', - help='fit under mpi4py.MPI.COMM_WORLD') - return p.parse_args() - - -def _load_cylinder_wake(data_path): - d = read_data(data_path) - dt = float(np.ravel(d['dt'])[0]) - u = np.asarray(d['u'], dtype=np.float64) - return u, dt - - -def _fit(args, comm): - u, dt = _load_cylinder_wake(args.data) - nt, n1, n2 = u.shape - is_root = comm is None or comm.rank == 0 - if is_root: - print(f'cylinder wake: nt={nt}, grid={n1}x{n2}, dt={dt}') - - data = u[..., np.newaxis] # single variable: streamwise velocity - params = dict( - n_dft=args.n_dft, - time_step=dt, - n_space_dims=2, - n_variables=1, - overlap=args.overlap, - regions=[1, 2], # sum- and difference-interactions - max_freq_idx=None, # the whole plane, as in the reference - solver='MengiOverton', - save_modes=False, # ~43k triads: modes would be ~8 GB - store_modes=False, - compute_energy_transfer=False, # T is not part of this figure - savedir=args.savedir, - ) - # refs/bmd/bmd.m:279-281 defaults to weight = ones(nx,1) ("uniform") when - # no weight is passed, and example1.m calls bmd(u) with none -- match that - # rather than a physically-motivated quadrature weight, since B (and so - # the colour scale) is linear in the weight - weights = utils_weights.uniform((n1, n2), n_vars=1, dV=1.0) - bmd = Standard(params=params, weights=weights, comm=comm).fit(data) - if is_root: - df = 1.0 / (args.n_dft * dt) - print(f'n_dft={args.n_dft} df={df:.6f} n_overlap={bmd.n_overlap} ' - f'n_blocks={bmd.n_blocks} n_triads={bmd.n_triads}') - return bmd - - -def _print_report(results): - triads = results.triads - values = np.abs(results.L[triads.f1_idx, triads.f2_idx]) - print('\nlabelled triads:') - for k, l in TRIADS: - i = triads.find(k, l) - print(f' (k,l,k+l) = ({k:3d},{l:3d},{k + l:3d}) ' - f'(f1,f2,f3) = ({triads.f1[i]:.4f}, {triads.f2[i]:.4f}, ' - f'{triads.f3[i]:.4f}) |lambda_1| = {values[i]:.4e}') - - print('\ntop 15 triads by |lambda_1| (k != 0 and l != 0):') - for row in top_triads(results, n=15, quantity='L', exclude_zero=True): - print(f" (k,l,k+l) = ({row['k']:3d},{row['l']:3d},{row['kl']:3d}) " - f"|lambda_1| = {row['value']:.4e}") - - -def _panel(ax, field, freq, df, xlim, ylim, cmap, vmin, vmax): - edges = np.append(freq - df / 2, freq[-1] + df / 2) - ax.pcolormesh(edges, edges, field.T, cmap=cmap, vmin=vmin, vmax=vmax, - shading='flat') - ax.set_aspect('equal') - ax.set_xlabel(r'$f_1$') - ax.set_ylabel(r'$f_2$') - ax.set_xlim(xlim) - ax.set_ylim(ylim) - - -def _plot(results, args): - freq, L = results.freq, results.L - df = freq[1] - freq[0] - triads = results.triads - field = np.ma.masked_invalid(np.log(np.abs(L))) - if args.clim is not None: - vmin, vmax = args.clim - else: - vmin, vmax = float(field.min()), float(field.max()) - cmap = _cylinder_colormap() - - fig, (axa, axb) = plt.subplots( - 1, 2, figsize=(7.6, 4.0), gridspec_kw=dict(width_ratios=[4, 3])) - - _panel(axa, field, freq, df, (0, freq[-1]), (-freq[-1], freq[-1] / 2), - cmap, vmin, vmax) - cax = axa.inset_axes([0.055, 0.026, 0.083, 0.251]) - fig.colorbar(ScalarMappable(Normalize(vmin, vmax), cmap), cax=cax, - ticks=[0, -10, -20], label=r'$\log(|\lambda_1|)$') - - _panel(axb, field, freq, df, (0, 0.8), (-0.8, 0.8), cmap, vmin, vmax) - - f0 = triads.f1[triads.find(12, 0)] # the shedding frequency, 12*df - axb.plot([0, 0.8], [f0, f0 - 0.8], 'k--', lw=0.8) - for k, l in TRIADS: - i = triads.find(k, l) - f1, f2 = triads.f1[i], triads.f2[i] - axb.plot(f1, f2, 'o', ms=7, mfc='none', mec='k', mew=1.0) - offset, ha = LABEL_LAYOUT[(k, l)] - axb.annotate(f'({k},{l})', (f1, f2), textcoords='offset points', - xytext=offset, fontsize=8, ha=ha, va='bottom') - - for ax, label in ((axa, '(a)'), (axb, '(b)')): - ax.text(-0.32, 1.08, label, transform=ax.transAxes, - fontsize=12, fontweight='bold', va='bottom') - - fig.tight_layout() - os.makedirs(os.path.dirname(args.out), exist_ok=True) - fig.savefig(args.out, dpi=args.dpi, bbox_inches='tight') - plt.close(fig) - print(f'\n[vmin, vmax] = [{vmin:.3f}, {vmax:.3f}]') - print(f'wrote {args.out}') - - -def main(): - args = _parse_args() - _check_prereqs(args.data) - - comm = None - if args.mpi: - from mpi4py import MPI - comm = MPI.COMM_WORLD - - if args.reuse: - results = load_results(args.savedir) - else: - _fit(args, comm) - if comm is not None and comm.rank != 0: - return - results = load_results(args.savedir) - - _print_report(results) - _plot(results, args) - - -if __name__ == '__main__': - main() diff --git a/docs/complex-conjugation-audit.md b/docs/complex-conjugation-audit.md new file mode 100644 index 0000000..6cbe62c --- /dev/null +++ b/docs/complex-conjugation-audit.md @@ -0,0 +1,98 @@ +# Complex-conjugation audit + +## Finding + +`normalize_data=True` is incorrect when PyBMD is given complex-valued time +domain data. In `Base._compute_blocks`, the variance estimate is calculated +as + +```python +np.sum((q_blk - np.mean(q_blk, axis=0))**2, axis=0) / den +``` + +([`pybmd/bmd/base.py:575`](../pybmd/bmd/base.py#L575)). For a complex random +variable, variance must use the Hermitian square, + +\[ + \operatorname{var}(q) = \frac{1}{N-1}\sum_t |q_t-\bar q|^2, +\] + +not the algebraic square \((q_t-\bar q)^2\). The latter is phase dependent +and can be complex or cancel to zero. The subsequent threshold and square +root then standardize by the wrong quantity. + +### Minimal reproducer + +```python +z = np.array([1, 1j, -1, -1j], dtype=complex) +d = z - z.mean() + +np.sum(d**2) / 3 # 0j: current PyBMD calculation +np.sum(np.abs(d)**2) / 3 # 4/3: correct variance +``` + +The current code replaces the zero result with one and leaves this signal with +variance `4/3`, rather than normalizing it to one. Other complex signals can +produce a complex scale factor, arbitrarily rotating and rescaling the data +before the FFT. This can change `B`, its numerical-radius optimizer result, +and the modes. Normal real-valued input is unaffected because `x**2` and +`abs(x)**2` agree for real `x`. + +## Comparison with the MATLAB reference + +The BMD matrix construction itself is *not* the fault: + +| Operation | PyBMD | `refs/bmd` | Result | +| --- | --- | --- | --- | +| Cross-spectral matrix | `q_sum.conj().T @ (...)` | `Q_hat_f3' * (...)` | Both use a conjugate transpose. | +| Mode norm | `np.vdot(psi, psi * w)` | `Psi' * (Psi .* weight)` | Both are Hermitian weighted inner products. | +| Transfer term | `np.vdot(psi_sum, psi_prod)` | `(...)' * (...)` | Both conjugate the sum-interaction factor. | + +MATLAB's apostrophe is conjugate transpose for complex arrays (as opposed to +`.'`, its non-conjugating transpose). Thus `refs/bmd/bmd.m:185,205,218` and +`refs/bmd/cbmd.m:161,181,196` use the right convention. PyBMD mirrors it in +`Base._triad_loop` (`base.py:646,663`) and `normalize_mode` +(`pybmd/bmd/utils.py:374`). Replacing any of these with a plain transpose +would be a separate, serious error, but that error is not present in the +audited BMD/CBMD matrix products. + +The reference has no `normalize_data` option, so it has no direct equivalent +of this defect. Consequently this is a PyBMD extension bug, not a PyBMD vs. +MATLAB porting mismatch. + +## Recommended fix and regression coverage + +Replace the variance line with a Hermitian magnitude square: + +```python +centered = q_blk - np.mean(q_blk, axis=0) +q_var = np.sum(np.abs(centered)**2, axis=0) / den +``` + +Add a regression test with the four-point reproducer above and assert that, +after normalization, each non-degenerate column has +`sum(abs(q_blk)**2) / (n_dft - 1) == 1`. A real-data regression should also +assert unchanged output, since the proposed expression is identical for real +inputs. + +## Verification performed + +- Inspected every conjugate-transpose/inner-product site in `pybmd/bmd` and + `refs/bmd/bmd.m` / `refs/bmd/cbmd.m`. +- Executed the reproducer: current variance `0j`; Hermitian variance + `1.3333333333333333`. +- Ran `pytest tests/optimizers -q`: **120 passed**. Those tests exercise the + complex numerical-radius matrices, but they do not cover complex input with + `normalize_data=True`. + +## Resolution + +Complex input never reached this code path in practice: `get_data_array` +cast the data to the requested float type, which silently discarded the +imaginary part (a `ComplexWarning` only). Both ends are now closed: + +- `pybmd/utils/io.py::get_data_array` raises `TypeError` on complex input, + since BMD's two-sided spectrum and sum/difference regions rely on the + conjugate symmetry of a real signal (`tests/test_io_rejects_complex.py`); +- the variance in `Base._compute_blocks` uses the Hermitian square + `np.abs(centered)**2`, which is identical for real data. diff --git a/docs/cylinder_bispectrum.md b/docs/cylinder_bispectrum.md deleted file mode 100644 index 1a8c23c..0000000 --- a/docs/cylinder_bispectrum.md +++ /dev/null @@ -1,81 +0,0 @@ -# Cylinder Wake Mode Bispectrum - -`docs/build_cylinder_figure.py` reproduces the cylinder-wake mode-bispectrum figure from Schmidt -(2020, *Nonlinear Dynamics*), `refs/figures/cylinder_bispectrum_sumdiff.pdf`, with PyBMD on the -`refs/bmd/wake_Re500.mat` dataset (the `refs/bmd` git submodule). - -```bash -git submodule update --init -MPLBACKEND=Agg python docs/build_cylinder_figure.py -``` - -![Cylinder wake mode bispectrum](figures/cylinder/cylinder_bispectrum_sumdiff.png) - -## How the reference parameters were recovered - -Nothing in the repo generates this figure -- there is no "sumdiff" script anywhere, and -`refs/bmd/example1.m` only draws an interactive, index-labelled version on the fly. The reference -PDF itself was rasterized at 400 dpi and measured directly (pcolor cell pitch, axis extents, marker -positions, colorbar ticks) to recover the parameters it was produced with: - -- panel (a) spans `xlim=[0, f(end)]`, `ylim=[-f(end), f(end)/2]` with `f(end) = 8.316`, i.e. - `df = 1/57.6 = 0.017361`; -- panel (b) zooms to `[0, 0.8] x [-0.8, 0.8]` and circles six triads at multiples of - `12 df = 0.2083`: `(12,12)`, `(12,0)`, `(24,12)`, `(24,24)`, `(36,12)`, `(36,24)`. - -Matching `df = 0.017361` requires `dt = 0.06`, `n_dft = 960` -- twice the time resolution of the -`wake_Re500.mat` shipped here (`dt = 0.12`, `nt = 1024`). The shipped file is that same dataset, -subsampled 2x in time over the same total duration (`nt * dt = 122.88` either way). Consequently: - -- **`n_dft = 480` at `dt = 0.12` reproduces the reference's exact frequency grid.** Verified with - `pybmd.bmd.utils.triad_indices`: all six labelled triads exist in `regions=[1, 2]` at exactly the - physical frequencies the reference marks them at, e.g. `(12,12,24)` -> `{0.2083, 0.2083, 0.4167}`. -- The dataset's own vortex-shedding frequency, measured from a zero-padded spectrum of `v`, is - `f0 = 0.2074` (harmonics at 0.4145, 0.6225, 0.8302) `= 11.94 df`, i.e. index **12** -- which is - why the reference labels the fundamental triad `(12,12,24)`. - -Two differences from the published panel are unavoidable given the shipped (subsampled) data, and -are not attempts to hide a bug: - -1. **Panel (a)'s extent is halved** -- `[0, 4.15] x [-4.15, 2.07]` instead of `[0, 8.32] x - [-8.32, 4.16]` -- because the Nyquist frequency halves when the sampling rate halves at fixed - `nt * dt`. Same picture, half the plane. -2. **3 blocks** at the reference's 50% overlap (`n_overlap = 240`), against presumably more in the - original full-rate dataset, so the non-resonant background sits higher and the lattice contrast - is a little weaker. - -The reference's colorbar (jet, with the dark-blue end faded to white) was reproduced by sampling -its pixels directly rather than guessed; see `_CBAR_STOPS` in `build_cylinder_figure.py`. - -### Spatial weight: match `bmd.m`'s own default, not a "better" one - -`B = Q3^H (Q1*Q2*w)/n_blocks` is linear in the spatial weight `w`, so `log|lambda_1|` (and the -whole colour scale) shifts by a constant depending on which weighting convention is used -- -independent of everything above. `refs/bmd/bmd.m:279-281` defaults to `weight = ones(nx,1)` -("uniform") when no weight is passed, and `refs/bmd/example1.m` calls `bmd(u)` with none, so the -published figure was made with a **uniform** weight, not a physically-motivated quadrature one. -Using `pybmd.utils.weights.trapz_2d` instead (a reasonable default for other PyBMD work) shifted -this figure's `[vmin, vmax]` to `[-29.99, -4.49]` against the reference's measured -`[-28.4, +0.37]`; switching to `pybmd.utils.weights.uniform((n1, n2), n_vars=1, dV=1.0)` (matching -`bmd.m`'s default) brings it to `[-25.9, -0.48]`. The residual ~1-in-log gap is consistent with the -reference likely bispectrum-ing `u` and `v` together (`n_variables=2` doubles the flattened -dimension and hence `|lambda_1|` roughly 2x) -- not reproduced here, since the figure only ever -labels a single field. - -## What to check when re-running this - -The script prints, for each of the six labelled triads, its resolved `(f1, f2, f3)` and `|lambda_1|`, -plus the top 15 triads overall (via `pybmd.bmd.postproc.top_triads`). The acceptance criterion is -physical, not pixel-exact: the labelled triads should sit among the strongest non-trivial -(`k != 0`, `l != 0`) entries, and the global maximum should fall on the shedding lattice (`k`, `l` -multiples of 12) -- exactly as observed: - -``` -top triads by |lambda_1| (k != 0 and l != 0): - (24,-12,12) 6.17e-01 - (12, 12,24) 6.11e-01 - (12,-12, 0) 2.29e-01 - ... -``` - -`(24,-12,12)` is the region-2 mirror of `(12,12,24)` and is expected to be comparably strong. diff --git a/docs/figures/cylinder/cylinder_bispectrum_sumdiff.png b/docs/figures/cylinder/cylinder_bispectrum_sumdiff.png deleted file mode 100644 index f433934..0000000 Binary files a/docs/figures/cylinder/cylinder_bispectrum_sumdiff.png and /dev/null differ diff --git a/docs/figures/hypothesis/hypothesis_harmonics_row.png b/docs/figures/hypothesis/hypothesis_harmonics_row.png deleted file mode 100644 index dc622d5..0000000 Binary files a/docs/figures/hypothesis/hypothesis_harmonics_row.png and /dev/null differ diff --git a/docs/figures/hypothesis/hypothesis_noise.png b/docs/figures/hypothesis/hypothesis_noise.png deleted file mode 100644 index a8acc5f..0000000 Binary files a/docs/figures/hypothesis/hypothesis_noise.png and /dev/null differ diff --git a/examples/example1_cylinder.py b/examples/example1_cylinder.py index 8025e70..76fe38e 100644 --- a/examples/example1_cylinder.py +++ b/examples/example1_cylinder.py @@ -14,7 +14,8 @@ sys.path.append(os.path.join(os.path.dirname(os.path.realpath(__file__)), '..')) from pybmd.bmd.standard import Standard -from pybmd.bmd.postproc import plot_mode_bispectrum, plot_triad_modes +from pybmd.bmd.postproc import (plot_mode_bispectrum, plot_triad_modes, + top_triads) import pybmd.utils.weights as utils_weights from examples.data import load_cylinder_wake @@ -42,21 +43,16 @@ def main(save_dir='example1_out'): weights = utils_weights.trapz_2d(x[:, 0], y[0, :], n_vars=1) bmd = Standard(params=params, weights=weights).fit(data) - ## the strongest triad, found through the triad map rather than by index - triads = bmd.triads - vals = np.abs(bmd.L[triads.f1_idx, triads.f2_idx]) - i_peak = int(np.argmax(vals[triads.k != 0])) - i_peak = int(np.flatnonzero(triads.k != 0)[i_peak]) - k, l = int(triads.k[i_peak]), int(triads.l[i_peak]) + ## the strongest triad with k, l both non-zero + peak = top_triads(bmd, n=1)[0] + k, l = int(peak['k']), int(peak['l']) print(f'strongest triad: (k,l,k+l) = ({k},{l},{k + l}), ' - f'|lambda_1| = {vals[i_peak]:.4e}') + f"|lambda_1| = {peak['value']:.4e}") - plot_mode_bispectrum( - bmd.L, bmd.freq, mark=[(triads.f1[i_peak], triads.f2[i_peak])], - path=save_dir, filename='bispectrum.png') - plot_triad_modes( - bmd.get_modes_at_triad(i_peak), k, l, x1=x[:, 0], x2=y[0, :], - path=save_dir, filename='modes.png') + plot_mode_bispectrum(bmd.L, bmd.freq, mark=[(peak['f1'], peak['f2'])], + path=save_dir, filename='bispectrum.png') + plot_triad_modes(bmd.get_modes_at_freqs(k, l), k, l, x1=x[:, 0], x2=y[0, :], + path=save_dir, filename='modes.png') print(f'figures written to {os.path.abspath(save_dir)}') return bmd diff --git a/examples/example4_hypothesis_testing.py b/examples/example4_hypothesis_testing.py new file mode 100644 index 0000000..83f4b64 --- /dev/null +++ b/examples/example4_hypothesis_testing.py @@ -0,0 +1,192 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +''' +Example 4: the "hypothesis testing" surrogate data of Schmidt (2020, Figs. 4 +and 5): travelling waves whose frequencies do or do not form a triad, with +and without noise. Renders the figures into ``example4_out/`` for visual +comparison with the paper. + + MPLBACKEND=Agg python examples/example4_hypothesis_testing.py + +``tests/test_hypothesis.py`` asserts the scientific content on the same data. +''' +import os +import sys + +import matplotlib.pyplot as plt +import numpy as np + +CFD = os.path.dirname(os.path.realpath(__file__)) +sys.path.append(os.path.join(CFD, '..')) + +from pybmd.bmd.standard import Standard +import pybmd.utils.weights as utils_weights + +# the paper's frequencies, moved to the nearest bins of n_dft=128 +NONRES = dict(name='nonres', freqs=(0.046875, 0.203125, 0.3515625)) # (0.05, 0.2, 0.35) +TRIAD = dict(name='triad', freqs=(0.046875, 0.203125, 0.25)) # (0.05, 0.2, 0.25) +QUARTET = dict(name='quartet', freqs=(0.046875, 0.1484375, 0.25, 0.453125)) # (0.05, 0.15, 0.25, 0.45) +NOISE = dict(name='noise', freqs=TRIAD['freqs'], snr=1.0) + + +def surrogate_waves(freqs, nt=1280, nx=100, dt=1.0, seed=0, snr=None): + ''' + ``q(x,t) = sum_j A_j cos(k_j x - 2 pi f_j t + theta0)``, unit amplitudes, + wavenumbers drawn from ``U[0, 5]`` on ``x in [0, 2 pi)`` with 100 points -- + exactly the paper's surrogate-data recipe. + + The paper adds a random phase offset per *realization*; a single + continuous time series segmented into 10 blocks of ``n_dft=128`` already + supplies that, since each block sees a different phase through ``t``, so + there is no need to simulate repeated realizations explicitly. ``snr=1`` + reproduces the paper's noise test: Gaussian noise scaled so its variance + equals the signal's. + ''' + rng = np.random.default_rng(seed) + x = np.linspace(0, 2 * np.pi, nx, endpoint=False) + t = np.arange(nt) * dt + k = rng.uniform(0, 5, size=len(freqs)) + q = np.zeros((nt, nx)) + for kj, fj in zip(k, freqs): + q += np.cos(kj * x[None, :] - 2 * np.pi * fj * t[:, None]) + if snr is not None: + q = q + rng.standard_normal(q.shape) * np.sqrt(q.var() / snr) + return q[..., np.newaxis], x, k + + +def fit_case(name, freqs, snr=None, save_dir='example4_out', + store_modes=False, **overrides): + '''BMD of one surrogate case with the paper's settings: n_dft=128, no + overlap, Hann window, sum interactions only; ``overrides`` go into + ``params``. Returns ``(bmd, x, k)``.''' + q, x, k = surrogate_waves(freqs, seed=0, snr=snr) + params = dict( + n_dft=128, time_step=1.0, n_space_dims=1, n_variables=1, overlap=0, + window='hann', regions=[1], solver='MengiOverton', save_modes=False, + store_modes=store_modes, savedir=os.path.join(save_dir, name)) + params.update(overrides) + w = utils_weights.uniform((x.size,), n_vars=1, dV=x[1] - x[0]) + return Standard(params=params, weights=w).fit(q), x, k + + +def _block_dft(q_x0, n_dft): + '''Hann-windowed DFT of the non-overlapping blocks of a single-point time + series, normalized like BMD's: ``(n_blocks, n_dft)``, column ``i`` being + integer frequency ``i`` (not fftshifted).''' + win = np.hanning(n_dft + 1)[:-1] + n_blocks = len(q_x0) // n_dft + blocks = (q_x0 - q_x0.mean())[:n_blocks * n_dft].reshape(n_blocks, n_dft) + return np.fft.fft(win * blocks, axis=1) / win.mean() / n_dft + + +def amplitude_spectrum(q_x0, n_dft, dt): + '''``A(f) = 2|mean_blocks q_hat(f)|``, computed independently of BMD with + the same window and blocking, as the paper's panel (a) does.''' + q_hat = _block_dft(q_x0, n_dft) + return np.fft.fftfreq(n_dft, dt), 2 * np.abs(q_hat).mean(axis=0) + + +def classical_bispectrum(q_x0, n_dft, m): + ''' + Classical (biased) bispectrum estimator of a single-point time series on + the integer-frequency grid ``0 <= j <= i < m``, block-averaged with the + same window and blocking as BMD -- the quantity the paper compares the + mode bispectrum against in its noise test. + ''' + q_hat = _block_dft(q_x0, n_dft) + B = np.full((m, m), np.nan) + for i in range(m): + for j in range(i + 1): + if i + j < n_dft // 2: + B[i, j] = np.abs(np.mean( + q_hat[:, i] * q_hat[:, j] * np.conj(q_hat[:, i + j]))) + return B + + +def _save(fig, path): + fig.savefig(path, dpi=150) + plt.close(fig) + print(f'wrote {path}') + + +def _surface(ax, t, vals, zmax, zlabel=r'$|\lambda_1|$'): + # plot_trisurf colours by the *data* range, not by set_zlim, so a panel + # that is flat relative to the z-axis would otherwise be painted with the + # full colormap and read as structured; pin vmin/vmax to the z-limits so + # colour and height agree, as MATLAB's fixed caxis does in the paper + ax.plot_trisurf(t.f1, t.f2, vals, cmap='viridis', linewidth=0.1, + vmin=0, vmax=zmax) + ax.set_zlim(0, zmax) + ax.set_xlabel('$f_1$') + ax.set_ylabel('$f_2$') + ax.set_zlabel(zlabel) + ax.set_xticks([0, 0.2, 0.4]) + ax.set_yticks([0, 0.1, 0.2]) + + +def _amplitude_panel(ax, q): + freq, A = amplitude_spectrum(q[:, 0, 0], 128, 1.0) + pos = freq >= 0 + ax.plot(freq[pos], A[pos], 'k') + ax.set_xlim(0, 0.5) + ax.set_xlabel('$f$') + + +def main(save_dir='example4_out'): + os.makedirs(save_dir, exist_ok=True) + + # -- figure 1: 3 rows (nonres / triad / quartet) x 2 columns ------------- + titles = { + 'nonres': r'$f_1 \pm f_2 \pm f_3 \neq 0$ (no triad)', + 'triad': r'$f_1 + f_2 = f_3$ (triad)', + 'quartet': (r'$f_1+f_2+f_3=f_4$, $f_k\pm f_l\pm f_m\neq 0$' + '\n(quartet, no triad)'), + } + colors = ['tab:blue', 'tab:red', 'tab:green', 'tab:purple'] + fig = plt.figure(figsize=(9, 12)) + for row, case in enumerate((NONRES, TRIAD, QUARTET)): + bmd, _, _ = fit_case(save_dir=save_dir, **case) + q, _, _ = surrogate_waves(case['freqs'], seed=0) + + ax_a = fig.add_subplot(3, 2, 2 * row + 1) + _amplitude_panel(ax_a, q) + for j, f in enumerate(case['freqs']): + ax_a.axvline(f, color=colors[j], lw=1) + ax_a.set_ylim(0, 1.05) + ax_a.set_ylabel('$A$') + ax_a.set_title(titles[case['name']], fontsize=9) + + t = bmd.triads + _surface(fig.add_subplot(3, 2, 2 * row + 2, projection='3d'), t, + np.abs(bmd.L[t.f1_idx, t.f2_idx]), 0.05) + fig.tight_layout() + _save(fig, os.path.join(save_dir, 'hypothesis_harmonics_row.png')) + + # -- figure 2: unit-SNR noise; amplitude, classical and mode bispectra --- + bmd, _, _ = fit_case(save_dir=save_dir, **NOISE) + q, _, _ = surrogate_waves(NOISE['freqs'], seed=0, snr=NOISE['snr']) + t = bmd.triads + B = classical_bispectrum(q[:, 0, 0], 128, 64) + + fig = plt.figure(figsize=(15, 4.5)) + ax0 = fig.add_subplot(1, 3, 1) + _amplitude_panel(ax0, q) + ax0.set_title('(a) amplitude spectrum') + + ax1 = fig.add_subplot(1, 3, 2, projection='3d') + _surface(ax1, t, B[t.k, t.l], 0.25, zlabel='$|B|$') + ax1.set_title('(b) classical bispectrum') + + ax2 = fig.add_subplot(1, 3, 3, projection='3d') + # the z-limit comes from the data (peak ~0.057): a fixed one either + # flattens the peak or leaves the panel mostly empty + vals = np.abs(bmd.L[t.f1_idx, t.f2_idx]) + _surface(ax2, t, vals, float(np.nanmax(vals)) * 1.05) + ax2.set_title('(c) mode bispectrum') + + fig.tight_layout() + _save(fig, os.path.join(save_dir, 'hypothesis_noise.png')) + + +if __name__ == '__main__': + main() diff --git a/examples/example5_cylinder_paper.md b/examples/example5_cylinder_paper.md new file mode 100644 index 0000000..5d0be4d --- /dev/null +++ b/examples/example5_cylinder_paper.md @@ -0,0 +1,135 @@ +# Cylinder Wake: Mode Bispectrum and Bispectral Modes + +`examples/example5_cylinder_paper.py` reproduces the cylinder-wake results of Schmidt (2020, +*Nonlinear Dynamics*) with PyBMD on the `refs/bmd/wake_Re500.mat` dataset (the `refs/bmd` git +submodule): the mode bispectrum of Fig. 7 (`refs/figures/cylinder_bispectrum_sumdiff.pdf`) and the +bispectral modes of Figs. 8 and 9 (see [below](#spatial-modes-figs-8-and-9)). + +```bash +git submodule update --init +MPLBACKEND=Agg python examples/example5_cylinder_paper.py # ~2-4.5 min, ~0.9 GB RAM +``` + +It writes `example5_out/cylinder_bispectrum_sumdiff.png` and `example5_out/modes_k{k}_l{l}.png` in +the working directory. Both use the same BMD settings with `q = [u, v]`, as in the paper. + +## How the reference parameters were recovered + +Nothing in the repo generates this figure -- there is no "sumdiff" script anywhere, and +`refs/bmd/example1.m` only draws an interactive, index-labelled version on the fly. The reference +PDF itself was rasterized at 400 dpi and measured directly (pcolor cell pitch, axis extents, marker +positions, colorbar ticks) to recover the parameters it was produced with: + +- panel (a) spans `xlim=[0, f(end)]`, `ylim=[-f(end), f(end)/2]` with `f(end) = 8.316`, i.e. + `df = 1/57.6 = 0.017361`; +- panel (b) zooms to `[0, 0.8] x [-0.8, 0.8]` and circles six triads at multiples of + `12 df = 0.2083`: `(12,12)`, `(12,0)`, `(24,12)`, `(24,24)`, `(36,12)`, `(36,24)`. + +Matching `df = 0.017361` requires `dt = 0.06`, `n_dft = 960` -- twice the time resolution of the +`wake_Re500.mat` shipped here (`dt = 0.12`, `nt = 1024`). The shipped file is that same dataset, +subsampled 2x in time over the same total duration (`nt * dt = 122.88` either way). Consequently: + +- **`n_dft = 480` at `dt = 0.12` reproduces the reference's exact frequency grid.** Verified with + `pybmd.bmd.utils.triad_indices`: all six labelled triads exist in `regions=[1, 2]` at exactly the + physical frequencies the reference marks them at, e.g. `(12,12,24)` -> `{0.2083, 0.2083, 0.4167}`. +- The dataset's own vortex-shedding frequency, measured from a zero-padded spectrum of `v`, is + `f0 = 0.2074` (harmonics at 0.4145, 0.6225, 0.8302) `= 11.94 df`, i.e. index **12** -- which is + why the reference labels the fundamental triad `(12,12,24)`. + +Two differences from the published panel are unavoidable given the shipped (subsampled) data, and +are not attempts to hide a bug: + +1. **Panel (a)'s extent is halved** -- `[0, 4.15] x [-4.15, 2.07]` instead of `[0, 8.32] x + [-8.32, 4.16]` -- because the Nyquist frequency halves when the sampling rate halves at fixed + `nt * dt`. Same picture, half the plane. +2. **3 blocks** at the reference's 50% overlap (`n_overlap = 240`), against presumably more in the + original full-rate dataset, so the non-resonant background sits higher and the lattice contrast + is a little weaker. + +The reference's colorbar (jet, with the dark-blue end faded to white) was reproduced by sampling +its pixels directly rather than guessed; see `CMAP` in `example5_cylinder_paper.py`. + +### Spatial weight: match `bmd.m`'s own default, not a "better" one + +`B = Q3^H (Q1*Q2*w)/n_blocks` is linear in the spatial weight `w`, so `log|lambda_1|` (and the +whole colour scale) shifts by a constant depending on which weighting convention is used -- +independent of everything above. `refs/bmd/bmd.m:279-281` defaults to `weight = ones(nx,1)` +("uniform") when no weight is passed, and `refs/bmd/example1.m` calls `bmd(u)` with none, so the +published figure was made with a **uniform** weight, not a physically-motivated quadrature one. +Using `pybmd.utils.weights.trapz_2d` instead (a reasonable default for other PyBMD work) shifted +this figure's `[vmin, vmax]` to `[-29.99, -4.49]` against the reference's measured +`[-28.4, +0.37]`; switching to `pybmd.utils.weights.uniform((n1, n2), n_vars=1, dV=1.0)` (matching +`bmd.m`'s default) brings it to `[-25.9, -0.48]` for `u` alone. + +The paper analyses `q = [u, v]`, so the script now uses both variables. That gives `[-25.55, -0.69]`, +which moves `vmax` slightly further from the reference rather than closing the gap. An earlier +version of this note guessed that the remaining difference of about 1 in `log|lambda_1|` came +from fitting `u` alone; this measurement rules that out. The remaining difference is unexplained. +The different number of blocks (item 2 above) is one candidate that has not been tested. + +## What to check when re-running this + +The script prints the top 15 triads overall (via `pybmd.bmd.postproc.top_triads`). The acceptance +criterion is physical, not pixel-exact: the labelled triads should sit among the strongest +non-trivial (`k != 0`, `l != 0`) entries, and the global maximum should fall on the shedding +lattice (`k`, `l` multiples of 12). A run with `q = [u, v]` gave: + +``` +top 15 triads (k != 0 and l != 0): + ( 12,-12, 0) |lambda_1| = 4.8380e-01 + ( 12, 12, 24) |lambda_1| = 4.7334e-01 + ( 24,-12, 12) |lambda_1| = 4.3104e-01 + ... +``` + +The paper places the global maximum at `(12,12,24)`. Here the difference self-interaction +`(12,-12,0)`, the `{f0, -f0, 0}` triad that drives the mean-flow deformation, is 2% above it. With +`u` alone the order was `(24,-12,12)` at 0.617 and then `(12,12,24)` at 0.611. In both cases the +leading triads lie on the shedding lattice and are within a few percent of each other. + +## Spatial modes (Figs. 8 and 9) + +The script does not redraw the paper's layout for these figures. It writes one `plot_triad_modes` +figure per triad circled in Fig. 7b, `example5_out/modes_k{k}_l{l}.png`, each showing the u and v +components of `phi_{k+l}`, `phi_{k o l}` and the interaction map `|phi_{k o l} phi_{k+l}|`. Fig. 8 +is the `phi_{k+l}` u panel of the six figures; Fig. 9 is the whole of `modes_k12_l12.png`. + +Holding the modes of all ~43k triads of the Fig. 7 fit would take ~15 GB. The modes therefore come +from a second fit with the same settings, restricted to `regions=[1]` and `max_freq_idx=36` (703 +triads). Each triad is solved independently, so the restriction leaves the results of the remaining +triads unchanged. The script checks this: `|lambda_1|` of the six labelled triads agrees between +the two fits to a relative difference of 0. + +The script prints the properties that Figs. 8 and 9 show. A run gave: + +``` + (k, l) n |lambda_1| sym(u) sym(v) lambda_x n*lambda_x + (12, 0) 1 5.0091e-01 0.001 1.000 3.998 3.998 + (12,12) 2 4.7334e-01 0.998 0.001 1.972 3.943 + (24,12) 3 9.3953e-02 0.001 0.999 1.371 4.112 + (24,24) 4 4.1311e-03 0.998 0.002 1.021 4.083 + (36,12) 4 2.0278e-02 0.996 0.005 1.023 4.091 + (36,24) 5 2.5207e-03 0.003 0.999 0.834 4.172 +Fig. 9: interaction map psi_{12,12} = |phi_{12o12} phi_{12+12}| + u: max 2.0877e-04 at (x, y) = (5.65, 0.65); 3.4% of the domain above half its max + v: max 8.7707e-04 at (x, y) = (12.54, 0.39); 11.6% of the domain above half its max +``` + +- `n` is the harmonic of `f0` at which `phi_{k+l}` oscillates. `sym` is the share of the mode's + energy that is even in `y`. For every triad, u is antisymmetric at odd `n` and symmetric at even + `n`, and v has the opposite parity. This is the symmetry of the vortex-shedding harmonics that + Fig. 8 shows. +- `lambda_x` is the peak streamwise wavelength of `phi_{k+l}`, measured on v. The product + `n*lambda_x` stays within 4.0 +/- 0.2 along the cascade, so each interaction shortens the + wavelength in proportion to the frequency. This matches the paper's statement that each + interaction yields new streamwise wavenumber components. +- `(12,12,24)` is the strongest triad with `k, l != 0` here too. +- For Fig. 9, the peak of the v interaction map is 4.2 times that of the u map, and the v map + covers 3.4 times as much of the domain above half its maximum. The paper reports both: "the + transverse component furthermore attains a larger maximum value than the streamwise component + and is less spatially confined". + +One claim of Fig. 9 cannot be checked on this dataset. The paper finds the interaction strongest +"in the wake region just downstream of the cylinder", but `wake_Re500.mat` covers only +`x = 2.7`-`14.9`, so the region near the cylinder is not in the data. The figures therefore show no +cylinder, and the location of the maximum above is the maximum within that window only. diff --git a/examples/example5_cylinder_paper.py b/examples/example5_cylinder_paper.py new file mode 100644 index 0000000..6e1ae11 --- /dev/null +++ b/examples/example5_cylinder_paper.py @@ -0,0 +1,204 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +''' +Example 5: the cylinder-wake results of Schmidt (2020, *Nonlinear Dynamics*) +on the full ``refs/bmd/wake_Re500.mat`` (``git submodule update --init``): +the mode bispectrum of Fig. 7, ``refs/figures/cylinder_bispectrum_sumdiff.pdf``, +and the bispectral modes of the triad cascade of Figs. 8 and 9. + + MPLBACKEND=Agg python examples/example5_cylinder_paper.py # ~4.5 min + +Writes to ``example5_out/``: + +- ``cylinder_bispectrum_sumdiff.png``, Fig. 7 in the paper's layout; +- ``modes_k*_l*.png``, one ``plot_triad_modes`` figure per triad circled in + Fig. 7b: the u and v components of the bispectral mode phi_{k+l}, the + cross-frequency field phi_{k o l} and the interaction map + |phi_{k o l} phi_{k+l}|. Fig. 8 is the phi_{k+l} row of all six; Fig. 9 is + the whole of ``modes_k12_l12.png``. + +Both use the same BMD settings, q = [u, v] as in the paper; the modes come +from a second fit restricted to the cascade (each triad is solved +independently, so the restriction changes nothing -- checked on lambda_1). +See example5_cylinder_paper.md for how the parameters were recovered. +''' +import os +import sys + +import matplotlib.pyplot as plt +from matplotlib.cm import ScalarMappable +from matplotlib.colors import LinearSegmentedColormap, Normalize +import numpy as np + +CFD = os.path.dirname(os.path.realpath(__file__)) +sys.path.append(os.path.join(CFD, '..')) + +from pybmd.bmd.standard import Standard +from pybmd.bmd.postproc import plot_triad_modes, top_triads +import pybmd.utils.weights as utils_weights +from pybmd.utils.io import read_data + +DATA_PATH = os.path.join(CFD, '..', 'refs', 'bmd', 'wake_Re500.mat') + +# the six triads circled in the reference's panel (b) -- the cascade of +# Fig. 8 -- with each label's offset (points) and alignment chosen so +# neighbouring labels don't collide +TRIADS = { + (12, 12): ((6, 2), 'left'), (12, 0): ((8, -2), 'left'), + (24, 12): ((-4, 8), 'right'), (24, 24): ((-4, 8), 'right'), + (36, 12): ((4, 8), 'left'), (36, 24): ((4, 8), 'left'), +} +F0_IDX = 12 # the shedding frequency f0, as an index + +# sampled from the reference's own colorbar: jet, with its dark-blue end +# faded to white +CMAP = LinearSegmentedColormap.from_list('cylinder_jet', ( + (0.000, (1.000, 1.000, 1.000)), (0.074, (0.835, 0.835, 1.000)), + (0.152, (0.482, 0.482, 1.000)), (0.230, (0.075, 0.075, 1.000)), + (0.307, (0.000, 0.345, 1.000)), (0.381, (0.004, 0.733, 1.000)), + (0.459, (0.051, 1.000, 0.953)), (0.537, (0.306, 1.000, 0.694)), + (0.615, (0.694, 1.000, 0.306)), (0.689, (1.000, 1.000, 0.000)), + (0.767, (1.000, 0.545, 0.000)), (0.844, (1.000, 0.090, 0.000)), + (0.922, (0.776, 0.000, 0.000)), (1.000, (0.502, 0.000, 0.000)), +)) + + +def _panel(ax, field, freq, xlim, ylim, vmin, vmax): + df = freq[1] - freq[0] + edges = np.append(freq - df / 2, freq[-1] + df / 2) + ax.pcolormesh(edges, edges, field.T, cmap=CMAP, vmin=vmin, vmax=vmax) + ax.set_aspect('equal') + ax.set_xlabel(r'$f_1$') + ax.set_ylabel(r'$f_2$') + ax.set_xlim(xlim) + ax.set_ylim(ylim) + + +def symmetric_fraction(field): + '''Share of ``sum |field|^2`` in the part even in y (axis 1), for a grid + symmetric about y = 0: 1 for a symmetric field, 0 for an antisymmetric one.''' + even = 0.5 * (field + field[:, ::-1]) + return float(np.sum(np.abs(even)**2) / np.sum(np.abs(field)**2)) + + +def streamwise_wavelength(field, dx, n_pad=4096): + '''Wavelength of the peak of the y-summed streamwise power spectrum of a + complex mode component, from a zero-padded DFT along x (axis 0).''' + spec = np.sum(np.abs(np.fft.fft(field, n=n_pad, axis=0))**2, axis=1) + kx = 2 * np.pi * np.fft.fftfreq(n_pad, dx) + spec[kx == 0] = 0 + return float(2 * np.pi / abs(kx[np.argmax(spec)])) + + +def _lambda1(bmd, k, l): + '''``|lambda_1|`` of the triad ``(k, l, k+l)`` of a fitted decomposition.''' + i = bmd.find_triad(k, l) + return abs(bmd.L[bmd.triads.f1_idx[i], bmd.triads.f2_idx[i]]) + + +def plot_bispectrum(bmd, path): + '''Fig. 7: (a) the sum and difference regions, (b) the low-frequency + magnification with the cascade circled.''' + triads, freq = bmd.triads, bmd.freq + field = np.ma.masked_invalid(np.log(np.abs(bmd.L))) + vmin, vmax = float(field.min()), float(field.max()) + fig, (axa, axb) = plt.subplots( + 1, 2, figsize=(7.6, 4.0), gridspec_kw=dict(width_ratios=[4, 3])) + + _panel(axa, field, freq, (0, freq[-1]), (-freq[-1], freq[-1] / 2), + vmin, vmax) + cax = axa.inset_axes([0.055, 0.026, 0.083, 0.251]) + fig.colorbar(ScalarMappable(Normalize(vmin, vmax), CMAP), cax=cax, + ticks=[0, -10, -20], label=r'$\log(|\lambda_1|)$') + + _panel(axb, field, freq, (0, 0.8), (-0.8, 0.8), vmin, vmax) + f0 = F0_IDX * (freq[1] - freq[0]) + axb.plot([0, 0.8], [f0, f0 - 0.8], 'k--', lw=0.8) + for (k, l), (offset, ha) in TRIADS.items(): + i = triads.find(k, l) + f1, f2 = triads.f1[i], triads.f2[i] + axb.plot(f1, f2, 'o', ms=7, mfc='none', mec='k', mew=1.0) + axb.annotate(f'({k},{l})', (f1, f2), textcoords='offset points', + xytext=offset, fontsize=8, ha=ha, va='bottom') + + for ax, label in ((axa, '(a)'), (axb, '(b)')): + ax.text(-0.32, 1.08, label, transform=ax.transAxes, + fontsize=12, fontweight='bold', va='bottom') + + fig.tight_layout() + fig.savefig(path, dpi=200, bbox_inches='tight') + plt.close(fig) + print(f'[vmin, vmax] = [{vmin:.3f}, {vmax:.3f}]') + print(f'wrote {path}') + + +def main(save_dir='example5_out'): + d = read_data(DATA_PATH) + dt = float(np.ravel(d['dt'])[0]) + x1, x2 = np.asarray(d['x'])[:, 0], np.asarray(d['y'])[0, :] + data = np.stack([d['u'], d['v']], axis=-1).astype(np.float64) + n1, n2, nv = data.shape[1:] + + params = dict( + n_dft=480, # the reference's df = 1/57.6 + time_step=dt, + n_space_dims=2, + n_variables=nv, # q = [u, v], as in the paper + overlap=50, + regions=[1, 2], # sum- and difference-interactions + solver='MengiOverton', + save_modes=False, + store_modes=False, # ~43k triads: modes would be ~15 GB + compute_energy_transfer=False, + savedir=os.path.join(save_dir, 'bispectrum'), + ) + # bmd.m defaults to a uniform weight and example1.m passes none; B is + # linear in the weight, so this sets the colour scale + weights = utils_weights.uniform((n1, n2), n_vars=nv, dV=1.0) + + # -- Fig. 7: mode bispectrum over the sum and difference regions --------- + bmd = Standard(params=params, weights=weights).fit(data) + print('top 15 triads (k != 0 and l != 0):') + for row in top_triads(bmd, n=15): + print(f" ({row['k']:3d},{row['l']:3d},{row['kl']:3d}) " + f"|lambda_1| = {row['value']:.4e}") + plot_bispectrum(bmd, os.path.join(save_dir, + 'cylinder_bispectrum_sumdiff.png')) + + # -- Figs. 8 and 9: modes, from the same settings on fewer triads -------- + params_modes = dict( + params, + regions=[1], # every triad of Fig. 8 is a sum-interaction + max_freq_idx=36, # covers the cascade, 703 triads + store_modes=True, # ~250 MB in memory + savedir=os.path.join(save_dir, 'modes'), + ) + bmd_modes = Standard(params=params_modes, weights=weights).fit(data) + + dx = x1[1] - x1[0] + print('Fig. 8: phi_{k+l}; sym = share of energy even in y ' + '(u: 0 for odd n, 1 for even n; v the opposite)') + print(' (k, l) n |lambda_1| sym(u) sym(v) lambda_x n*lambda_x') + max_rel = 0.0 + for k, l in TRIADS: + lam1 = _lambda1(bmd, k, l) + max_rel = max(max_rel, abs(_lambda1(bmd_modes, k, l) - lam1) / lam1) + modes = bmd_modes.get_modes_at_freqs(k, l) + phi = modes[0] + n = (k + l) // F0_IDX + lam = streamwise_wavelength(phi[..., 1], dx) + print(f' ({k:2d},{l:2d}) {n:2d} {lam1:.4e} ' + f'{symmetric_fraction(phi[..., 0]):.3f} ' + f'{symmetric_fraction(phi[..., 1]):.3f} ' + f'{lam:7.3f} {n * lam:7.3f}') + plot_triad_modes( + modes, k, l, x1=x1, x2=x2, vars_idx=(0, 1), figsize=(11, 9), + xlabel='$x$', ylabel='$y$', path=save_dir, + filename=f'modes_k{k}_l{l}.png') + print(f'max relative difference of |lambda_1| between the two fits: ' + f'{max_rel:.1e}') + return bmd, bmd_modes + + +if __name__ == '__main__': + main() diff --git a/pybmd/bmd/CLAUDE.md b/pybmd/bmd/CLAUDE.md index 1eb8fbe..c05e3aa 100644 --- a/pybmd/bmd/CLAUDE.md +++ b/pybmd/bmd/CLAUDE.md @@ -6,22 +6,28 @@ commands and the project overview. ## The Base/Standard/Cross split -`Base` owns everything: params, weights, mean, DFT blocking, the triad map, the solve loop, MPI, -and storage. **Subclasses override only three things**, so the algorithm lives in exactly one place: +`Base` owns everything, including `fit()` itself: params, weights, mean, DFT blocking, the triad +map, the solve loop, MPI, and storage. **Subclasses override only the per-triad matrices and a few +shape hooks**, so the algorithm lives in exactly one place: | | `Standard` (BMD) | `Cross` (CBMD) | | --- | --- | --- | | `_triad_matrices(q_hat, i)` | returns `(Q3, Q1*Q2, weights)` | stacks `n_state` blocks; sums the `q*r` terms | | `_constituent_matrices(q_hat, i)` | returns `(Q1, Q2)` (only with `constituent_modes`) | not supported — rejected at construction | -| `_compute_qhat` block shape | `(nx*nv,)` | `(nx, nv)` so `q_hat[f][:, v]` is contiguous | -| `define_weights` | `(*xshape, nv)` | overridden: `xshape`, **no variable axis** | +| `_block_shape()` (one `q_hat` row) | `(nx*nv,)` | `(nx, nv)` so `q_hat[f][:, v]` is contiguous | +| `_mode_shape` (property) | `(*xshape, nv)` | `(*xshape, n_state)` | +| `_expected_weights_shape()` | `(*xshape, nv)` | `xshape`, **no variable axis** | +| `_post_initialize()` | no-op | tiles the weights over the states, prints `state_idx`/`qr_idx` | +| `_unflatten_modes(psi)` | plain reshape | state-slowest unflatten, see below | +| `_label` | `'BMD'` | `'CBMD'` (timing print only) | Everything after `B = Q_sum^H (Q_prod * w) / n_blocks` is shared in `Base._triad_loop`. ## Data flow -`fit()` → `_initialize` (dims, `n_blocks`, weights, mean, `Triads`, savedir — clearing stale -`modes/triad_idx_*.npy` from an earlier run into the same directory — size guard) → +`Base.fit()` → `_initialize` (dims, `n_blocks`, weights, mean, `Triads`, savedir — clearing stale +`modes/triad_idx_*.npy` from an earlier run into the same directory — size guard, then +`_post_initialize`) → `_compute_qhat` → `_triad_loop` → `_store_and_save` (arrays first, `params_modes.yaml` last, so a YAML failure cannot lose results). @@ -37,8 +43,8 @@ some triad actually references. With `max_freq_idx` set that is a small fraction the states along the flat axis with the state **slowest** (`flat = j*nx + p`, matching `cbmd.m`'s `repmat`), so a C-order reshape straight into `(*xshape, n_state)` scrambles every `n_state > 1` mode. `_unflatten_modes` is the single place a flat mode becomes a field — - `Cross` overrides it to unflatten as `(n_state, *xshape)` and move the state axis last — - and `tests/test_cbmd.py::test_two_identical_states_give_identical_mode_slices` guards it. + `Cross` overrides it to unflatten as `(n_state, *xshape)` and move the state axis last; + `tests/test_octave_reference.py`'s CBMD tiers compare the modes against `cbmd.m`. The other hazards are the weights and a user `mean`: both are therefore checked against the full `(*xshape, nv)` shape and a bare flat vector (or the reference's variable-first layout) is **rejected**. @@ -77,7 +83,7 @@ some triad actually references. With `max_freq_idx` set that is a small fraction Each fixes a silent wrong answer; all three are covered by regression tests. Measured end-to-end on the 169 triads of the full cylinder-wake dataset (`regions=[1,2]`, `max_freq_idx=12`), run *directly under Octave* against `refs/bmd/bmd.m` itself (see -[`docs/octave_cross_validation.md`](../../docs/octave_cross_validation.md)): `MengiOverton` matches a +[`tests/octave/octave_cross_validation.md`](../../tests/octave/octave_cross_validation.md)): `MengiOverton` matches a brute-force scan of the numerical radius to ~5e-8, the genuine `refs/bmd.m` is off by >1% on 52 of the 169 triads (>10% on 29), always an *under*-estimate, since `B = Q3^H (Q1∘Q2∘w)/n_blocks` is tiny (median `‖B‖₁ ~ 5.2e-6` there). These figures were originally measured against a Python @@ -108,7 +114,7 @@ the original), otherwise its absolute `|w − w_old| ≤ tol` stopping test fire the tiny matrices BMD produces. Confirmed live under Octave, for both `bmd.m` and `cbmd.m` (see -[`docs/octave_cross_validation.md`](../../docs/octave_cross_validation.md)): the reference's actually +[`tests/octave/octave_cross_validation.md`](../../tests/octave/octave_cross_validation.md)): the reference's actually *reachable* solvers are `'MengiOverton'` and `'HeWatson'`. `'simpleIteration'` passes the option validator but the inner `switch` has no matching case (`case {'simpleit'}` is what's there instead) and errors with `'Unknown solver.'`; `'eig'` fails the same way; `'simpleit'` itself fails the @@ -130,7 +136,7 @@ first pass. PyBMD has no use for a solver that needs an unseeded random start ve `simpleIteration` already reproduces the underlying algorithm deterministically via `default_start` — confirmed live on the paper's own hypothesis-test triad case (`tests/test_hypothesis.py`'s surrogate-data recipe, run through both implementations; see -[`docs/octave_cross_validation.md`](../../docs/octave_cross_validation.md)): `simpleIteration` agrees +[`tests/octave/octave_cross_validation.md`](../../tests/octave/octave_cross_validation.md)): `simpleIteration` agrees with `MengiOverton` everywhere there (max relative deviation 4.4e-4, 0/780 triads above 1%), while `refs/bmd.m`'s `HeWatson` disagrees with both on up to 753/780 triads on the flat, non-resonant case — random-start non-convergence on a featureless @@ -155,7 +161,7 @@ without noise (`‖B‖₁ ~ 3.5e-9`) — reproducing a branch decision taken ex boundary is inherently unstable across LAPACK builds and language boundaries, and is not a defect to chase further. Validated live against Octave in `tests/test_octave_reference.py::test_tier_c_matlab_compat_reproduces_reference`; see -[`docs/octave_cross_validation.md`](../../docs/octave_cross_validation.md) for the figures. +[`tests/octave/octave_cross_validation.md`](../../tests/octave/octave_cross_validation.md) for the figures. ## Conventions that bite diff --git a/pybmd/bmd/base.py b/pybmd/bmd/base.py index 031772d..ebe8817 100644 --- a/pybmd/bmd/base.py +++ b/pybmd/bmd/base.py @@ -1,9 +1,9 @@ ''' -Base module for the BMD: - - The `Base.fit` method must be implemented in inherited classes +Base module for the BMD: parameters, weights, mean, DFT blocking, the triad +loop and storage. :class:`~pybmd.bmd.standard.Standard` and +:class:`~pybmd.bmd.cross.Cross` only supply the per-triad matrices and the +shape hooks below. ''' -from __future__ import division - import glob import os import time @@ -57,6 +57,8 @@ class Base(): computed from the data. Default is None. ''' + _label = 'BMD' # name used in the timing print + def __init__(self, params, weights=None, comm=None, mean=None): ##--- required self._n_dft = params['n_dft'] @@ -143,17 +145,43 @@ def __init__(self, params, weights=None, comm=None, mean=None): self._window = self._set_dtype(self._window) self._resolve_overlap() - # -------------------------------------------------------------------------- - # to be implemented by inherited classes - # -------------------------------------------------------------------------- - - def fit(self, data_list, *args, **kwargs): + def fit(self, data_list): ''' - Fit the data using BMD. + Fit the data: initialize, DFT every block, solve every triad, save. + + :param data_list: data matrix of shape ``(nt, *xshape, n_variables)``, + or path(s) to it. - :param list data_list: data matrix for which to compute the BMD. + :return: the fitted object. ''' - raise NotImplementedError # pragma: no cover + start0 = time.time() + + start = time.time() + self._initialize(data_list) + self._pr0(f'Time to initialize: {time.time() - start} s.') + + start = time.time() + q_hat = self._compute_qhat() + self._pr0(f'Time to compute DFT: {time.time() - start} s.') + del self.data + utils_par.barrier(self._comm) + + start = time.time() + self._triad_loop(q_hat) + del q_hat + self._pr0(f'------------------------------------') + self._pr0(f'Time to compute {self._label}: {time.time() - start} s.') + + self._store_and_save() + self._pr0(f' ') + self._pr0(f'Results saved in folder {self._savedir_sim}') + self._pr0(f'Total time: {time.time() - start0} s.') + utils_par.barrier(self._comm) + return self + + # -------------------------------------------------------------------------- + # hooks for inherited classes + # -------------------------------------------------------------------------- def _triad_matrices(self, q_hat, i_triad): ''' @@ -182,10 +210,19 @@ def _expected_weights_shape(self): '''Shape a user weight array must have: spatial shape plus variables.''' return tuple(self._xshape) + (self._nv,) - def _mode_elements(self): - '''Number of entries of one mode, i.e. the length of the flat axis - that :meth:`_triad_matrices` builds.''' - return self._nxv + @property + def _mode_shape(self): + '''Shape of one mode: the flat axis :meth:`_triad_matrices` builds, + unflattened.''' + return (*self._xshape, self._nv) + + def _block_shape(self): + '''Shape of one frequency row of a block, as stored in ``q_hat``.''' + return (self._nxv,) + + def _post_initialize(self): + '''Subclass set-up that needs the flattened weights; runs last in + :meth:`_initialize`.''' def _unflatten_modes(self, psi): ''' @@ -402,8 +439,11 @@ def _initialize(self, data_list): self._n_blocks = num // den # test feasibility - if (self._n_dft < 4) or (self._n_blocks < 2): - raise ValueError('Spectral estimation parameters not meaningful.') + if self._n_blocks < 2: + raise ValueError( + f'Spectral estimation parameters not meaningful: nt={self._nt}, ' + f'n_dft={self._n_dft}, n_overlap={self._n_overlap} give ' + f'{self._n_blocks} block(s), at least 2 are needed.') ## define and check weights self.define_weights() @@ -462,7 +502,7 @@ def _initialize(self, data_list): * self._n_blocks * self._complex(1).nbytes * B2GB) self._modes_size_gb = (self._n_mode_comp * self.n_triads - * self._mode_elements() + * int(np.prod(self._mode_shape)) * self._complex(1).nbytes * B2GB) if ((self._save_modes or self._store_modes) and self._modes_size_gb > self._max_modes_gb): @@ -476,6 +516,7 @@ def _initialize(self, data_list): self._print_parameters() self._pr0(f'------------------------------------') + self._post_initialize() def define_weights(self): '''Define and check weights.''' @@ -572,7 +613,8 @@ def _compute_blocks(self, i_blk): # standardize every point and variable to unit variance within the # block, i.e. divide by the standard deviation den = self._n_dft - 1 - q_var = np.sum((q_blk - np.mean(q_blk, axis=0))**2, axis=0) / den + centered = q_blk - np.mean(q_blk, axis=0) + q_var = np.sum(np.abs(centered)**2, axis=0) / den q_var[q_var < 4 * np.finfo(q_blk.dtype).eps] = 1 q_blk = q_blk / np.sqrt(q_var) @@ -581,23 +623,22 @@ def _compute_blocks(self, i_blk): q_blk_hat = (self._win_weight / self._n_dft) * np.fft.fft(q_blk, axis=0) return np.fft.fftshift(q_blk_hat, axes=0), offset - def _compute_qhat(self, block_shape): + def _compute_qhat(self): ''' Fourier realizations for every frequency row any triad refers to. Only the rows in ``triads.freq_needed`` are retained; for a bispectrum restricted by ``max_freq_idx`` that is a small fraction of ``n_dft``. - :param tuple block_shape: shape of one frequency row of a block. - :return: mapping from frequency row to its ``(*block_shape, n_blocks)`` - array of realizations. + array of realizations, ``block_shape`` being :meth:`_block_shape`. :rtype: dict ''' self._pr0(f' ') self._pr0(f'Calculating temporal DFT') self._pr0(f'------------------------------------') + block_shape = self._block_shape() needed = self._triads.freq_needed q_hat = {int(f): np.empty((*block_shape, self._n_blocks), dtype=self._complex) for f in needed} diff --git a/pybmd/bmd/cross.py b/pybmd/bmd/cross.py index 05c89e9..b9d952c 100644 --- a/pybmd/bmd/cross.py +++ b/pybmd/bmd/cross.py @@ -1,10 +1,7 @@ '''Derived module from base.py for cross-bispectral mode decomposition.''' -import time - import numpy as np from pybmd.bmd.base import Base -import pybmd.utils.parallel as utils_par class Cross(Base): @@ -37,6 +34,8 @@ class Cross(Base): Nonlinear Dynamics, 2020. DOI 10.1007/s11071-020-06037-z ''' + _label = 'CBMD' + def __init__(self, params, weights=None, comm=None, mean=None): super().__init__(params, weights=weights, comm=comm, mean=mean) self._state_idx = np.atleast_1d( @@ -93,8 +92,20 @@ def _expected_weights_shape(self): ''' return tuple(self._xshape) - def _mode_elements(self): - return self.n_state * self._nx + @property + def _mode_shape(self): + return (*self._xshape, self.n_state) + + def _block_shape(self): + # (nx, nv), so that q_hat[f][:, v] is a contiguous slice + return (self._nx, self._nv) + + def _post_initialize(self): + # the same spatial weight applies to every state; tile the whole + # spatial vector n_state times (a per-element repeat would scramble it) + self._weights_tiled = np.tile(self._weights, (self.n_state, 1)) + self._pr0(f'State indices : {self._state_idx.tolist()}') + self._pr0(f'q*r indices : {self._qr_idx.tolist()}') def _unflatten_modes(self, psi): ''' @@ -108,48 +119,6 @@ def _unflatten_modes(self, psi): psi = psi.reshape((psi.shape[0], self.n_state, *self._xshape)) return np.moveaxis(psi, 1, -1) - def fit(self, data_list): - ''' - Class-specific method to fit the data matrix using the CBMD algorithm. - - :param data_list: data matrix of shape ``(nt, *xshape, n_variables)``, - or path(s) to it. - - :return: the fitted object. - :rtype: Cross - ''' - start0 = time.time() - - start = time.time() - self._initialize(data_list) - self._mode_shape = (*self._xshape, self.n_state) - # the same spatial weight applies to every state; tile the whole - # spatial vector n_state times (a per-element repeat would scramble it) - self._weights_tiled = np.tile(self._weights, (self.n_state, 1)) - assert self._weights_tiled.shape == (self.n_state * self._nx, 1) - self._pr0(f'State indices : {self._state_idx.tolist()}') - self._pr0(f'q*r indices : {self._qr_idx.tolist()}') - self._pr0(f'Time to initialize: {time.time() - start} s.') - - start = time.time() - q_hat = self._compute_qhat(block_shape=(self._nx, self._nv)) - self._pr0(f'Time to compute DFT: {time.time() - start} s.') - del self.data - utils_par.barrier(self._comm) - - start = time.time() - self._triad_loop(q_hat) - del q_hat - self._pr0(f'------------------------------------') - self._pr0(f'Time to compute CBMD: {time.time() - start} s.') - - self._store_and_save() - self._pr0(f' ') - self._pr0(f'Results saved in folder {self._savedir_sim}') - self._pr0(f'Total time: {time.time() - start0} s.') - utils_par.barrier(self._comm) - return self - def _triad_matrices(self, q_hat, i_triad): ''' Assemble the state realizations and the quadratic term for one triad, diff --git a/pybmd/bmd/optimizers.py b/pybmd/bmd/optimizers.py index 3577672..793bd8b 100644 --- a/pybmd/bmd/optimizers.py +++ b/pybmd/bmd/optimizers.py @@ -259,7 +259,7 @@ def mengi_overton(A, tol=1e-8, n_it_max=500, matlab_compat=False): weights), every genuine level-set crossing is rejected and the search returns a local value at ``theta=0`` -- always an *under*-estimate, confirmed live under Octave against the real ``bmd.m``/``cbmd.m`` (see - ``docs/octave_cross_validation.md``): 52/169 triads off by >1% (29 by + ``tests/octave/octave_cross_validation.md``): 52/169 triads off by >1% (29 by >10%) on the full cylinder-wake fixture. ``matlab_compat=True`` reproduces that under-estimate to ~4e-6 relative when ``B`` is reasonably well scaled (``||B||_1 >~ 1e-6``); as ``||B||_1`` falls @@ -267,7 +267,7 @@ def mengi_overton(A, tol=1e-8, n_it_max=500, matlab_compat=False): noise-free cases of the paper's hypothesis test) the branch decisions it is reproducing sit exactly at the tolerance boundary, so agreement degrades and is not a defect to chase further -- see - ``docs/octave_cross_validation.md`` for the measured figures. + ``tests/octave/octave_cross_validation.md`` for the measured figures. :return: the value ``w = z^H A z`` and the maximiser ``z``. :rtype: tuple(complex, numpy.ndarray) diff --git a/pybmd/bmd/postproc.py b/pybmd/bmd/postproc.py index 04dfef5..6c1be6a 100644 --- a/pybmd/bmd/postproc.py +++ b/pybmd/bmd/postproc.py @@ -155,8 +155,8 @@ def top_triads(results, n=10, quantity='L', exclude_zero=True): ''' Return the strongest triads in a loaded result. - :param BMDResults or str results: loaded results, or a directory accepted - by :func:`load_results`. + :param results: a :class:`BMDResults`, a fitted ``Standard``/``Cross``, + or a results directory accepted by :func:`load_results`. :param int n: number of triads to return. :param str quantity: ``'L'`` for mode bispectrum or ``'T'`` for energy transfer magnitude. @@ -166,7 +166,7 @@ def top_triads(results, n=10, quantity='L', exclude_zero=True): frequencies, region, and value. :rtype: numpy.ndarray ''' - if not isinstance(results, BMDResults): + if isinstance(results, (str, os.PathLike)): results = load_results(results) quantity = quantity.upper() diff --git a/pybmd/bmd/standard.py b/pybmd/bmd/standard.py index fce058d..fe55c29 100644 --- a/pybmd/bmd/standard.py +++ b/pybmd/bmd/standard.py @@ -1,8 +1,5 @@ '''Derived module from base.py for standard BMD.''' -import time - from pybmd.bmd.base import Base -import pybmd.utils.parallel as utils_par class Standard(Base): @@ -21,42 +18,6 @@ class Standard(Base): Nonlinear Dynamics, 2020. DOI 10.1007/s11071-020-06037-z ''' - def fit(self, data_list): - ''' - Class-specific method to fit the data matrix using the BMD algorithm. - - :param data_list: data matrix of shape ``(nt, *xshape, n_variables)``, - or path(s) to it. - - :return: the fitted object. - :rtype: Standard - ''' - start0 = time.time() - - start = time.time() - self._initialize(data_list) - self._mode_shape = (*self._xshape, self._nv) - self._pr0(f'Time to initialize: {time.time() - start} s.') - - start = time.time() - q_hat = self._compute_qhat(block_shape=(self._nxv,)) - self._pr0(f'Time to compute DFT: {time.time() - start} s.') - del self.data - utils_par.barrier(self._comm) - - start = time.time() - self._triad_loop(q_hat) - del q_hat - self._pr0(f'------------------------------------') - self._pr0(f'Time to compute BMD: {time.time() - start} s.') - - self._store_and_save() - self._pr0(f' ') - self._pr0(f'Results saved in folder {self._savedir_sim}') - self._pr0(f'Total time: {time.time() - start0} s.') - utils_par.barrier(self._comm) - return self - def _triad_matrices(self, q_hat, i_triad): ''' Assemble the realizations of the sum interaction and of the quadratic diff --git a/pybmd/utils/io.py b/pybmd/utils/io.py index 8da6f19..e076ad1 100644 --- a/pybmd/utils/io.py +++ b/pybmd/utils/io.py @@ -1,13 +1,8 @@ '''Module implementing I/O utils used across the library.''' -import argparse import os from os.path import splitext import numpy as np -import yaml - - -REQUIRED_KEYS = ['time_step', 'n_space_dims', 'n_variables', 'n_dft'] def read_data(data_file, format=None, comm=None): @@ -100,37 +95,6 @@ def _as_array(obj): return np.asarray(obj) -def read_config(parsed_file=None): - ''' - Parse a YAML config file with ``required:`` and ``optional:`` sections, - each a list of single-key mappings. - - :param str parsed_file: file to parse. Default is None, in which case the - path is read from the ``--config_file`` command-line argument. - - :return: the parameters read from the config file. - :rtype: dict - ''' - parser = argparse.ArgumentParser(description='Config file.') - parser.add_argument('--config_file', required=True, - help='Configuration file.') - if parsed_file: - args = parser.parse_args(['--config_file', parsed_file]) - else: - args = parser.parse_args() - - with open(args.config_file) as file: - l = yaml.load(file, Loader=yaml.FullLoader) - - params = _parse_yaml(l['required']) - found, missing = _check_keys(params, REQUIRED_KEYS) - if not found: - raise ValueError(f'config file is missing required keys: {missing}') - if 'optional' in l: - params = {**params, **_parse_yaml(l['optional'])} - return params - - def get_data_array(data_list, xdim, nv, dtype=np.float64): ''' Assemble the input into a single array of shape ``(nt, *xshape, nv)``. @@ -174,19 +138,10 @@ def get_data_array(data_list, xdim, nv, dtype=np.float64): raise ValueError( f'data has {data.shape[-1]} variables in its last axis but ' f'n_variables is {nv}.') + if np.iscomplexobj(data): + # casting to float below would silently drop the imaginary part + raise TypeError( + 'PyBMD expects real-valued data: the two-sided spectrum and the ' + 'sum/difference regions rely on the conjugate symmetry of a real ' + 'signal. Pass the real part explicitly if that is what you mean.') return np.ascontiguousarray(data, dtype=dtype) - - -def _parse_yaml(l): - params = dict() - for d in l: - k = list(d.keys())[0] - params[k] = d[k] - return params - - -def _check_keys(l, keys): - if isinstance(keys, str): - keys = [keys] - keys_not_found = [k for k in keys if k not in l.keys()] - return not keys_not_found, keys_not_found diff --git a/pybmd/utils/parallel.py b/pybmd/utils/parallel.py index 9318164..56ad61c 100644 --- a/pybmd/utils/parallel.py +++ b/pybmd/utils/parallel.py @@ -60,56 +60,25 @@ def allreduce(data, comm): return reduced -def allreduce_scalar(value, comm, op='sum'): +def distribute_indices(n, comm): ''' - Reduce a scalar across all ranks. - - :param value: local value. - :param MPI.Comm comm: parallel communicator, or None. - :param str op: one of 'sum', 'min', 'max'. Default is 'sum'. - - :return: the reduced value, identical on every rank. - ''' - if comm is None: - return value - MPI = _get_module_MPI(comm) - ops = {'sum': MPI.SUM, 'min': MPI.MIN, 'max': MPI.MAX} - return comm.allreduce(value, op=ops[op]) - - -def _blockdist(n, size, rank): - '''Contiguous block distribution of ``n`` items; returns (count, start).''' - q, r = divmod(n, size) - count = q + (1 if r > rank else 0) - start = rank * q + min(rank, r) - return (count, start) if rank < size else (0, 0) - - -def distribute_indices(n, comm, mode='round_robin'): - ''' - Split ``range(n)`` across ranks. + Split ``range(n)`` across ranks, round-robin. :param int n: number of items, here the number of triads. :param MPI.Comm comm: parallel communicator, or None. - :param str mode: 'round_robin' (default) or 'block'. :return: the indices owned by this rank. :rtype: numpy.ndarray .. note:: - The default is round-robin rather than contiguous blocks because the - cost of a triad varies systematically across the ``f1``-``f2`` plane: - the numerical-radius solve takes more iterations where the spectrum of + Round-robin rather than contiguous blocks because the cost of a triad + varies systematically across the ``f1``-``f2`` plane: the + numerical-radius solve takes more iterations where the spectrum of ``B`` is clustered, which happens in bands. A contiguous split would hand one rank an entire band; interleaving balances the load with an imbalance of at most one triad. ''' if comm is None: return np.arange(n) - if mode == 'round_robin': - return np.arange(comm.rank, n, comm.size) - if mode == 'block': - count, start = _blockdist(n, comm.size, comm.rank) - return np.arange(start, start + count) - raise ValueError(f"mode must be 'round_robin' or 'block'; got {mode!r}.") + return np.arange(comm.rank, n, comm.size) diff --git a/pybmd/utils/weights.py b/pybmd/utils/weights.py index 5fbd439..a195965 100644 --- a/pybmd/utils/weights.py +++ b/pybmd/utils/weights.py @@ -83,14 +83,6 @@ def trapz_3d(x1, x2, x3, n_vars=1): return {'weights_name': 'trapz_3d', 'weights': dV} -def custom(**kwargs): - ''' - Customized weights, to be implemented by the user if required. The returned - array must have the shape documented at the top of this module. - ''' - pass - - def apply_normalization(data, weights, n_vars, method='variance', comm=None): ''' Normalize the weights variable-wise by the data variance. diff --git a/tests/CLAUDE.md b/tests/CLAUDE.md index 6a63816..88e06f1 100644 --- a/tests/CLAUDE.md +++ b/tests/CLAUDE.md @@ -3,30 +3,28 @@ Testing strategy detail. See the root [`CLAUDE.md`](../CLAUDE.md) for the pytest commands and [`pybmd/bmd/CLAUDE.md`](../pybmd/bmd/CLAUDE.md) for the solver deviations these tests guard. -The primary regression net does not depend on MATLAB and runs everywhere — `test_bmd_serial.py`, -`test_cbmd.py`, `test_bmd_mpi.py`, `test_io.py`: +Four layers, from solver to paper: -- `test_bispectrum_matches_closed_form`: an on-grid, boxcar-windowed, non-overlapping signal with - components at `k1`, `k2`, `k1+k2` carrying independent random phases per block, the third - locked to the sum of the first two. Every `B` is then rank one, so for *every* triad - `L = (a_k·a_l·a_{k+l}/8)·Σ w·conj(φ_{k+l})·φ_k·φ_l` (non-zero only for the 12 sum/difference - views of that one physical triad — `(5,3)`, `(8,-3)`, `(8,-5)` in regions 1, 2 — and exactly - zero elsewhere), `T` is the same sum without the weights, and the modes are the spatial - patterns themselves. Matches to ~1e-12 relative. `expected_bispectrum` in that file is the - oracle; extend it rather than hard-coding values. -- conjugate symmetry `L(-k,-l) = conj(L(k,l))`, `T(-k,-l) = T(k,l)` on random data, all 8 regions; -- exact triad counts: `(m+1)²` for regions {1,2}, `(m+1)(m+2)/2` for {1}, a brute-force count of - `|k+l| < Nyquist` for all eight (625 / 325 / 2401 / **12223** at `n_dft=128` — not 12288); -- CBMD reducing to BMD triad-for-triad when q=r=s, and two identical CBMD states giving identical - mode slices and `L = 2·L_BMD` — the guard for `Cross._unflatten_modes` (state-slowest flat axis); -- `test_bmd_mpi.py`: bit-identical `L`, `T`, `coeffs` and every mode file between `mpirun -n 1` - and `-n 2` (marker `mpi`, self-skips without `mpirun`/`mpi4py`); its helper also checks - `allreduce` on a big-endian buffer; -- `test_io.py`: the h5py (MATLAB v7.3) and scipy (v5) `.mat` paths return the same arrays. +- `tests/optimizers/` (one test function per file): the numerical-radius solvers against a + brute-force angular scan and against each other — signed `max_fov`, `_pow2_scale` (including + subnormals), scale equivariance, the `MengiOvertonMATLAB` under-estimate, iteration caps, + determinism. These run in a few seconds and need nothing but NumPy/SciPy. +- `test_hypothesis.py`: Schmidt (2020)'s hypothesis test on the surrogate data of + `examples/example4_hypothesis_testing.py` (whose `fit_case`/`surrogate_waves` it imports) — the + resonant triad is detected at the right bin and scale, non-resonant and quartet cases stay + flat, the peak survives unit-SNR noise, and the modes recover the travelling waves. Four fits + are memoized across the module. +- `test_bmd_mpi.py` (marker `mpi`, self-skips without `mpirun`/`mpi4py`): bit-identical `L`, `T`, + `coeffs` and every mode file between `mpirun -n 1` and `-n 2`, through `tests/mpi_fit.py`; its + helper also checks `allreduce` on a big-endian buffer. +- `test_octave_reference.py` (marker `slow`): the reference `refs/bmd/bmd.m`/`cbmd.m` run live + under Octave, in three tiers — A: `Q_hat` and every per-triad `B` from an instrumented copy, + isolating the DFT/blocking/weighting stage from the solver; B: `L`, `T` and the modes end to + end; C: the solver deviations and `MengiOvertonMATLAB`'s bug-compatibility, measured. + `tests/octave/octave_ref.py` is the harness, `tests/octave/build_report.py` regenerates the + figures of [`octave/octave_cross_validation.md`](octave/octave_cross_validation.md). +- `test_io_rejects_complex.py`: complex data is refused rather than cast to its real part. -Every bug-fix regression in these files (weights mutation, CBMD mode scrambling, `normalize_data`, -stale mode files, mean layout, `_pow2_scale` subnormals, `simple_iteration` tolerance, the h5py -transpose, the signed `T` plot, ...) was confirmed to fail on the code before its fix. `tests/conftest.py` puts the checkout first on `sys.path`, so a non-editable install cannot shadow the source. @@ -34,17 +32,15 @@ Modes are defined only up to a unit-modulus phase — compare them with `|| / (‖a‖‖b‖) ≈ 1`, never elementwise. Octave is available on this machine and can run `refs/bmd/bmd.m`/`cbmd.m` directly, so the -Deviations numbers in [`pybmd/bmd/CLAUDE.md`](../pybmd/bmd/CLAUDE.md) are now cross-checked against +Deviations numbers in [`pybmd/bmd/CLAUDE.md`](../pybmd/bmd/CLAUDE.md) are cross-checked against the genuine MATLAB source rather than only a Python transcription of it. `refs/bmd` is a **git submodule** pointing at `olivertschmidt/bmd` — the reference's research/non-commercial license means it must stay a pointer rather than vendored code; run `git submodule update --init` to -populate it locally. `tests/test_octave_reference.py` (`pytest -m slow`, or `pytest -tests/test_octave_reference.py`) exercises this; it self-skips, not errors, when `octave-cli` is +populate it locally. `test_octave_reference.py` self-skips, not errors, when `octave-cli` is absent or the submodule hasn't been initialized. One caveat: `bmd.m`'s Hamming window `hammwin` is a file-local subfunction Octave cannot call from outside `bmd.m`, so `test_default_window_matches_reference` evaluates the *transcribed formula* under Octave — it -checks NumPy against Octave arithmetic on `bmd.m:310`'s expression, not the file itself. `.github/workflows/octave_reference.yml` runs it -in CI (checks out the submodule, `apt-get install`s Octave). See -[`docs/octave_cross_validation.md`](../docs/octave_cross_validation.md) for the method (three -comparison tiers, isolating the DFT/blocking/weighting stage from the solver) and the full measured -tables, with figures. +checks NumPy against Octave arithmetic on `bmd.m:310`'s expression, not the file itself. +`.github/workflows/octave_reference.yml` runs it in CI (checks out the submodule, `apt-get +install`s Octave). See [`octave/octave_cross_validation.md`](octave/octave_cross_validation.md) +for the method and the full measured tables, with figures. diff --git a/tests/data/input_bmd.yaml b/tests/data/input_bmd.yaml deleted file mode 100644 index 526a48c..0000000 --- a/tests/data/input_bmd.yaml +++ /dev/null @@ -1,17 +0,0 @@ -required: - - time_step : 1 - - n_space_dims: 2 - - n_variables : 1 - - n_dft : 32 - -optional: - - overlap : 50 - - mean_type : 'longtime' - - regions : [1, 2] - - max_freq_idx : 8 - - solver : 'MengiOverton' - - tol : 1.0e-6 - - n_it_max : 500 - - dtype : 'double' - - savedir : 'bmd_results' - - save_modes : True diff --git a/docs/build_octave_report.py b/tests/octave/build_report.py similarity index 82% rename from docs/build_octave_report.py rename to tests/octave/build_report.py index 838f4d6..b0c6454 100644 --- a/docs/build_octave_report.py +++ b/tests/octave/build_report.py @@ -1,13 +1,14 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- ''' -Regenerate the figures embedded in ``docs/octave_cross_validation.md``. +Regenerate the figures embedded in ``octave_cross_validation.md`` (written +to ``figures/`` next to this file). Requires ``octave-cli`` on PATH and the ``refs/bmd`` submodule populated (``git submodule update --init``); run from anywhere, paths are resolved relative to this file. - python docs/build_octave_report.py + python tests/octave/build_report.py ''' import os import shutil @@ -19,12 +20,11 @@ import numpy as np import scipy.io -DOCS_DIR = os.path.dirname(os.path.realpath(__file__)) -REPO_ROOT = os.path.realpath(os.path.join(DOCS_DIR, '..')) -FIG_DIR = os.path.join(DOCS_DIR, 'figures', 'octave') +OCTAVE_DIR = os.path.dirname(os.path.realpath(__file__)) +REPO_ROOT = os.path.realpath(os.path.join(OCTAVE_DIR, '..', '..')) +FIG_DIR = os.path.join(OCTAVE_DIR, 'figures') sys.path.insert(0, REPO_ROOT) -sys.path.insert(0, os.path.join(REPO_ROOT, 'tests', 'octave')) -sys.path.insert(0, os.path.join(REPO_ROOT, 'tests')) +sys.path.insert(0, OCTAVE_DIR) from pybmd.bmd.standard import Standard from pybmd.bmd.postproc import plot_mode_bispectrum @@ -32,9 +32,22 @@ import pybmd.utils.weights as utils_weights import octave_ref as oref -# the paper's surrogate-data recipe, reused rather than duplicated -- see -# test_hypothesis.py's own docstring for the recipe itself -from test_hypothesis import surrogate_waves, TRIAD +# the paper's surrogate-data recipe and BMD settings, reused rather than +# duplicated +from examples.example4_hypothesis_testing import (surrogate_waves, fit_case, + TRIAD) + + +def _save(fig, name): + path = os.path.join(FIG_DIR, name) + fig.savefig(path, dpi=140) + plt.close(fig) + print(f'wrote {path}') + + +def _rel(a, ref): + '''Element-wise relative deviation of ``a`` from ``ref``.''' + return np.abs(a - ref) / np.maximum(ref, 1e-300) def _check_prereqs(): @@ -51,7 +64,7 @@ def _check_prereqs(): def _full_dataset_run(pybmd_solver='MengiOverton'): ''' - PyBMD and reference L, at the config docs/octave_cross_validation.md + PyBMD and reference L, at the config octave_cross_validation.md cites. ``pybmd_solver`` selects PyBMD's own solver; the reference side always runs bmd.m's own MengiOverton -- the fixed point of comparison. ''' @@ -87,10 +100,7 @@ def fig_bispectrum_comparison(bmd, L_ref): fig.suptitle('Mode bispectrum $\\log|\\lambda_1|$ -- cylinder wake, ' 'regions={1,2}, max_freq_idx=12, 169 triads') fig.tight_layout() - path = os.path.join(FIG_DIR, 'bispectrum_comparison.png') - fig.savefig(path, dpi=140) - plt.close(fig) - print(f'wrote {path}') + _save(fig, 'bispectrum_comparison.png') def fig_deviation(bmd, L_ref): @@ -98,7 +108,7 @@ def fig_deviation(bmd, L_ref): t = bmd.triads vals_py = np.abs(bmd.L[t.f1_idx, t.f2_idx]) vals_ref = np.abs(L_ref[t.f1_idx, t.f2_idx]) - rel = np.abs(vals_ref - vals_py) / np.maximum(vals_py, 1e-300) + rel = _rel(vals_ref, vals_py) fig, ax = plt.subplots(figsize=(7, 6.5)) sc = ax.scatter(t.k, t.l, c=100 * rel, cmap='inferno_r', s=45, @@ -116,10 +126,7 @@ def fig_deviation(bmd, L_ref): ax.legend(loc='upper right', frameon=True, fontsize=9) fig.colorbar(sc, ax=ax, label='relative deviation, %') fig.tight_layout() - path = os.path.join(FIG_DIR, 'deviation_heatmap.png') - fig.savefig(path, dpi=140) - plt.close(fig) - print(f'wrote {path}') + _save(fig, 'deviation_heatmap.png') return rel @@ -165,7 +172,7 @@ def fig_three_way_solver_comparison(bmd, L_ref): (5, vals_ref, 'Reference deviation\nfrom PyBMD default, %'), (6, vals_compat, 'MengiOvertonMATLAB deviation\nfrom PyBMD default, %')): ax = fig.add_subplot(2, 3, ax_idx) - rel = np.abs(vals - vals_py) / np.maximum(vals_py, 1e-300) + rel = _rel(vals, vals_py) sc = ax.scatter(t.k, t.l, c=100 * rel, cmap='inferno_r', s=40, vmin=0, vmax=max(1.0, float(100 * rel.max())), edgecolors='none') @@ -178,14 +185,11 @@ def fig_three_way_solver_comparison(bmd, L_ref): fig.suptitle('Three-way solver comparison -- cylinder wake, regions={1,2}, ' 'max_freq_idx=12, 169 triads') fig.tight_layout() - path = os.path.join(FIG_DIR, 'three_way_solver_comparison.png') - fig.savefig(path, dpi=140) - plt.close(fig) - print(f'wrote {path}') + _save(fig, 'three_way_solver_comparison.png') - rel_ref = np.abs(vals_ref - vals_py) / np.maximum(vals_py, 1e-300) - rel_compat = np.abs(vals_compat - vals_py) / np.maximum(vals_py, 1e-300) - rel_compat_vs_ref = np.abs(vals_compat - vals_ref) / np.maximum(vals_ref, 1e-300) + rel_ref = _rel(vals_ref, vals_py) + rel_compat = _rel(vals_compat, vals_py) + rel_compat_vs_ref = _rel(vals_compat, vals_ref) print(f' bmd.m vs PyBMD default: max rel {rel_ref.max():.3e}; ' f'>1%: {int((rel_ref > 0.01).sum())}/{len(rel_ref)}; ' f'>10%: {int((rel_ref > 0.10).sum())}/{len(rel_ref)}') @@ -195,31 +199,26 @@ def fig_three_way_solver_comparison(bmd, L_ref): f'>1%: {int((rel_compat_vs_ref > 0.01).sum())}/{len(rel_compat_vs_ref)}') -def _hypothesis_run(freqs, snr, max_freq_idx=40, n_dft=128, seed=0): +def _hypothesis_run(freqs, snr, max_freq_idx=40): ''' PyBMD (three solvers) and the reference bmd.m (two solvers) on one - hypothesis-test surrogate case -- see test_hypothesis.py's - ``surrogate_waves`` for the recipe this reproduces exactly (n_dft=128, + hypothesis-test surrogate case, through example4's ``fit_case`` (n_dft=128, overlap=0, Hann window, regions=[1], 10 blocks). :return: ``(results, triads)``, where ``results`` maps - ``'pybmd_'`` to a fitted :class:`Standard` and - ``'bmd_'`` to the reference's raw ``L``. + ``'pybmd_'`` and ``'bmd_'`` to the respective ``L``. ''' - q, x, k = surrogate_waves(freqs, seed=seed, snr=snr) + q, x, k = surrogate_waves(freqs, seed=0, snr=snr) w = utils_weights.uniform((x.size,), n_vars=1, dV=x[1] - x[0]) results = {} for solver in ('MengiOverton', 'MengiOvertonMATLAB', 'simpleIteration'): - params = dict(n_dft=n_dft, time_step=1.0, n_space_dims=1, n_variables=1, - overlap=0, window='hann', regions=[1], - max_freq_idx=max_freq_idx, solver=solver, save_modes=False, - savedir=os.path.join(FIG_DIR, f'_scratch_hyp_{solver}')) - bmd = Standard(params=params, weights=w).fit(q) - results[f'pybmd_{solver}'] = bmd - shutil.rmtree(params['savedir'], ignore_errors=True) - - bmd = results['pybmd_MengiOverton'] + name = f'_scratch_hyp_{solver}' + bmd, _, _ = fit_case(name, freqs, snr=snr, save_dir=FIG_DIR, + solver=solver, max_freq_idx=max_freq_idx) + results[f'pybmd_{solver}'] = bmd.L + shutil.rmtree(os.path.join(FIG_DIR, name), ignore_errors=True) + # bmd.m's HeWatson draws an unseeded random start vector (refs/bmd/bmd.m # has no seeding hook this driver can reach), so its numbers -- unlike # every other figure in this script -- vary run to run; that variability @@ -251,7 +250,7 @@ def fig_hypothesis_pybmd_vs_matlab(): results, t = _hypothesis_run(TRIAD['freqs'], snr) panels = [ ('PyBMD MengiOverton', - np.abs(results['pybmd_MengiOverton'].L[t.f1_idx, t.f2_idx])), + np.abs(results['pybmd_MengiOverton'][t.f1_idx, t.f2_idx])), ('bmd.m MengiOverton', np.abs(results['bmd_MengiOverton'][t.f1_idx, t.f2_idx])), ('bmd.m HeWatson', @@ -269,7 +268,7 @@ def fig_hypothesis_pybmd_vs_matlab(): ax.set_title(f'{name} ({label})\npeak ({t.k[i]},{t.l[i]}) ' f'$|\\lambda_1|$={vals[i]:.5f}', fontsize=8) - py = np.abs(results['pybmd_MengiOverton'].L[t.f1_idx, t.f2_idx]) + py = np.abs(results['pybmd_MengiOverton'][t.f1_idx, t.f2_idx]) order = np.argsort(py) ax = fig.add_subplot(2, 4, row * 4 + 4) alternatives = [ @@ -279,9 +278,8 @@ def fig_hypothesis_pybmd_vs_matlab(): ('PyBMD MengiOvertonMATLAB', 'pybmd_MengiOvertonMATLAB', 'tab:blue'), ] for name, key, color in alternatives: - vals = (np.abs(results[key].L[t.f1_idx, t.f2_idx]) if key.startswith('pybmd_') - else np.abs(results[key][t.f1_idx, t.f2_idx])) - rel = np.abs(vals - py) / np.maximum(py, 1e-300) + vals = np.abs(results[key][t.f1_idx, t.f2_idx]) + rel = _rel(vals, py) ax.semilogy(np.arange(len(py)), np.maximum(rel[order], 1e-16), '.', ms=3, color=color, label=name) summary.append((label, name, float(rel.max()), @@ -300,10 +298,7 @@ def fig_hypothesis_pybmd_vs_matlab(): # a manual rect plus explicit spacing avoids the overlap tight_layout # alone leaves between them. fig.subplots_adjust(top=0.86, bottom=0.08, hspace=0.45, wspace=0.5) - path = os.path.join(FIG_DIR, 'hypothesis_pybmd_vs_matlab.png') - fig.savefig(path, dpi=140) - plt.close(fig) - print(f'wrote {path}') + _save(fig, 'hypothesis_pybmd_vs_matlab.png') for label, name, mx, n1, n in summary: print(f' [{label}] {name:24s} vs PyBMD MengiOverton: ' f'max rel {mx:.3e}; >1%: {n1}/{n}') @@ -343,11 +338,8 @@ def fig_scale_equivariance(): ax.set_aspect('equal') ax.legend(loc='upper left', frameon=True) fig.tight_layout() - path = os.path.join(FIG_DIR, 'scale_equivariance.png') - fig.savefig(path, dpi=140) - plt.close(fig) - print(f'wrote {path}') - rel = np.abs(b - a) / np.maximum(a, 1e-300) + _save(fig, 'scale_equivariance.png') + rel = _rel(b, a) print(f' max rel deviation from equivariance: {rel.max():.3f}; ' f'>1%: {int((rel > 0.01).sum())}/{rel.size}; ' f'>10%: {int((rel > 0.10).sum())}/{rel.size}') diff --git a/docs/figures/octave/bispectrum_comparison.png b/tests/octave/figures/bispectrum_comparison.png similarity index 100% rename from docs/figures/octave/bispectrum_comparison.png rename to tests/octave/figures/bispectrum_comparison.png diff --git a/docs/figures/octave/deviation_heatmap.png b/tests/octave/figures/deviation_heatmap.png similarity index 100% rename from docs/figures/octave/deviation_heatmap.png rename to tests/octave/figures/deviation_heatmap.png diff --git a/docs/figures/octave/hypothesis_pybmd_vs_matlab.png b/tests/octave/figures/hypothesis_pybmd_vs_matlab.png similarity index 100% rename from docs/figures/octave/hypothesis_pybmd_vs_matlab.png rename to tests/octave/figures/hypothesis_pybmd_vs_matlab.png diff --git a/docs/figures/octave/scale_equivariance.png b/tests/octave/figures/scale_equivariance.png similarity index 100% rename from docs/figures/octave/scale_equivariance.png rename to tests/octave/figures/scale_equivariance.png diff --git a/docs/figures/octave/three_way_solver_comparison.png b/tests/octave/figures/three_way_solver_comparison.png similarity index 100% rename from docs/figures/octave/three_way_solver_comparison.png rename to tests/octave/figures/three_way_solver_comparison.png diff --git a/docs/octave_cross_validation.md b/tests/octave/octave_cross_validation.md similarity index 95% rename from docs/octave_cross_validation.md rename to tests/octave/octave_cross_validation.md index 40e4759..d1ffebc 100644 --- a/docs/octave_cross_validation.md +++ b/tests/octave/octave_cross_validation.md @@ -61,20 +61,20 @@ Side by side on the full cylinder-wake fixture (`regions=[1,2]`, `max_freq_idx=1 PyBMD's mode bispectrum and the reference's look qualitatively the same but disagree exactly where Tier C predicts: -![Mode bispectrum: PyBMD vs. the reference](figures/octave/bispectrum_comparison.png) +![Mode bispectrum: PyBMD vs. the reference](figures/bispectrum_comparison.png) The per-triad relative deviation, mapped onto the `(k,l)` plane, never exceeds PyBMD and is concentrated where `|λ₁|` is smallest — an under-estimate confined to the weak triads, not a uniform mismatch: -![Per-triad deviation heatmap](figures/octave/deviation_heatmap.png) +![Per-triad deviation heatmap](figures/deviation_heatmap.png) A correct solver for the numerical radius is exactly scale-equivariant, `r(cA) = c·r(A)`. Running the unmodified reference on a random case and on a `1e-2` rescale of it shows the reference itself violating this by up to 32% on this random case — independent confirmation that the fault is in `bmd.m`'s solver, not in anything PyBMD does to the data before comparing: -![Reference scale-equivariance error](figures/octave/scale_equivariance.png) +![Reference scale-equivariance error](figures/scale_equivariance.png) ## Solver comparison and the MATLAB-compatible solver @@ -123,7 +123,7 @@ essentially the whole gap, and all three only reproduce `bmd.m` applied together | only `_pow2_scale` reverted | 2.810e-01 | 1/81 | | only the `max(w,1)` clamp reverted | 6.663e-01 | 7/81 | -![Three-way solver comparison](figures/octave/three_way_solver_comparison.png) +![Three-way solver comparison](figures/three_way_solver_comparison.png) The deviation maps in the (k,l) plane (bottom row) are visually near-identical between the reference and `MengiOvertonMATLAB` — the compat solver reproduces not just the aggregate counts @@ -154,7 +154,7 @@ surrogate-data recipe in `tests/test_hypothesis.py` (`n_dft=128`, `overlap=0`, H | triad, no noise | 3.48e-09 | (26, 6), 0.04052 | | triad, SNR = 1 | 2.15e-03 | (26, 6), 0.05659 | -![Hypothesis test, PyBMD vs. bmd.m](figures/octave/hypothesis_pybmd_vs_matlab.png) +![Hypothesis test, PyBMD vs. bmd.m](figures/hypothesis_pybmd_vs_matlab.png) Per-triad deviation from PyBMD's `MengiOverton`: @@ -183,10 +183,10 @@ Three conclusions: ## Regenerating Figures -The report figures in `docs/figures/octave/` can be regenerated with: +The report figures in `tests/octave/figures/` can be regenerated with: ```bash -python docs/build_octave_report.py +python tests/octave/build_report.py ``` That script requires Octave and the populated `refs/bmd` submodule, and takes a few minutes: it diff --git a/tests/octave/octave_ref.py b/tests/octave/octave_ref.py index 6e719e5..f584367 100644 --- a/tests/octave/octave_ref.py +++ b/tests/octave/octave_ref.py @@ -2,7 +2,7 @@ Cross-validation harness against the original MATLAB ``bmd.m``/``cbmd.m``, run unmodified (Tier B) or lightly instrumented (Tier A) under Octave. -See ``docs/octave_cross_validation.md`` for the fairness rules this module +See ``tests/octave/octave_cross_validation.md`` for the fairness rules this module exists to enforce (weight flatten order, single- vs double-precision input, window/overlap parity, the CBMD variable-axis position, ...), and ``tests/test_octave_reference.py`` for the tests that use it. diff --git a/tests/optimizers/test_matlab_compat_underestimates.py b/tests/optimizers/test_matlab_compat_underestimates.py index 734087d..4c61cde 100644 --- a/tests/optimizers/test_matlab_compat_underestimates.py +++ b/tests/optimizers/test_matlab_compat_underestimates.py @@ -31,7 +31,7 @@ def test_matlab_compat_underestimates_on_tiny_matrix(): w_compat, _ = mengi_overton(A, tol=1e-6, n_it_max=500, matlab_compat=True) # never exceeds the corrected solver, matching the one-sided invariant - # measured live against the real bmd.m (docs/octave_cross_validation.md) + # measured live against the real bmd.m (tests/octave/octave_cross_validation.md) assert abs(w_compat) <= abs(w_default) * (1 + 1e-8) # and, at this scale, is materially lower -- not merely numerically equal assert abs(w_compat) < abs(w_default) * 0.9 diff --git a/tests/test_hypothesis.py b/tests/test_hypothesis.py index 46f4fb1..29ae3a3 100644 --- a/tests/test_hypothesis.py +++ b/tests/test_hypothesis.py @@ -1,6 +1,11 @@ #!/usr/bin/env python3 # -*- coding: utf-8 -*- -'''Validation against the "hypothesis testing" surrogate data of Schmidt (2020).''' +'''Validation against the "hypothesis testing" surrogate data of Schmidt (2020). + +The surrogate data and spectral estimators live in +``examples/example4_hypothesis_testing.py``, which also renders the reference +figures. +''' import atexit import os import shutil @@ -14,53 +19,13 @@ CFD = os.path.dirname(CF) sys.path.append(os.path.join(CFD, '../')) -from pybmd.bmd.standard import Standard -import pybmd.utils.weights as utils_weights - -FIGURES_DIR = os.path.join(CFD, '..', 'docs', 'figures', 'hypothesis') +from examples.example4_hypothesis_testing import ( + NONRES, TRIAD, QUARTET, NOISE, surrogate_waves, fit_case, + amplitude_spectrum, classical_bispectrum) _TMPDIR = tempfile.mkdtemp(prefix='pybmd_test_bmd_') atexit.register(shutil.rmtree, _TMPDIR, ignore_errors=True) - -# --------------------------------------------------------------------------- -# surrogate data -- "Hypothesis testing" section of Schmidt (2020) -# --------------------------------------------------------------------------- - -def surrogate_waves(freqs, nt=1280, nx=100, dt=1.0, seed=0, snr=None): - ''' - ``q(x,t) = sum_j A_j cos(k_j x - 2 pi f_j t + theta0)``, unit amplitudes, - wavenumbers drawn from ``U[0, 5]`` on ``x in [0, 2 pi)`` with 100 points -- - exactly the paper's surrogate-data recipe. - - The paper adds a random phase offset per *realization*; a single - continuous time series segmented into 10 blocks of ``n_dft=128`` already - supplies that, since each block sees a different phase through ``t``, so - there is no need to simulate repeated realizations explicitly. ``snr=1`` - reproduces the paper's noise test: Gaussian noise scaled so its variance - equals the signal's. - ''' - rng = np.random.default_rng(seed) - x = np.linspace(0, 2 * np.pi, nx, endpoint=False) - t = np.arange(nt) * dt - k = rng.uniform(0, 5, size=len(freqs)) - q = np.zeros((nt, nx)) - for kj, fj in zip(k, freqs): - q += np.cos(kj * x[None, :] - 2 * np.pi * fj * t[:, None]) - if snr is not None: - q = q + rng.standard_normal(q.shape) * np.sqrt(q.var() / snr) - return q[..., np.newaxis], x, k - - -def _params(name, **kwargs): - params = dict( - n_dft=128, time_step=1.0, n_space_dims=1, n_variables=1, overlap=0, - window='hann', regions=[1], solver='MengiOverton', save_modes=False, - savedir=os.path.join(_TMPDIR, name)) - params.update(kwargs) - return params - - _CASES = {} @@ -68,56 +33,11 @@ def _case(name, freqs, **kwargs): '''Fit one paper case, memoized -- every test in this module shares the same 4 fits (~2-4 s each), rather than re-fitting per assertion.''' if name not in _CASES: - q, x, k = surrogate_waves(freqs, seed=0, **kwargs) - w = utils_weights.uniform((x.size,), n_vars=1, dV=x[1] - x[0]) - bmd = Standard(params=_params(name, store_modes=(name == 'triad')), - weights=w).fit(q) - _CASES[name] = (bmd, x, k) + _CASES[name] = fit_case(name, freqs, save_dir=_TMPDIR, + store_modes=(name == 'triad'), **kwargs) return _CASES[name] -NONRES = dict(name='nonres', freqs=(0.046875, 0.203125, 0.3515625)) # (0.05, 0.2, 0.35) -TRIAD = dict(name='triad', freqs=(0.046875, 0.203125, 0.25)) # (0.05, 0.2, 0.25) -QUARTET = dict(name='quartet', freqs=(0.046875, 0.1484375, 0.25, 0.453125)) # (0.05, 0.15, 0.25, 0.45) - - -def amplitude_spectrum(q_x0, n_dft, dt, window='hann'): - '''``A(f) = 2|mean_blocks q_hat(f)|``, computed independently of BMD with - the same window and blocking, as the paper's panel (a) does.''' - win = np.hanning(n_dft + 1)[:-1] - win_weight = 1.0 / win.mean() - n_blocks = len(q_x0) // n_dft - q_c = q_x0 - q_x0.mean() - blocks = np.stack([q_c[i * n_dft:(i + 1) * n_dft] for i in range(n_blocks)]) - q_hat = np.fft.fft(win[None, :] * blocks, axis=1) * win_weight / n_dft - freq = np.fft.fftfreq(n_dft, dt) - return freq, 2 * np.abs(q_hat).mean(axis=0) - - -def classical_bispectrum(q_x0, n_dft, dt, m, window='hann'): - ''' - Classical (biased) bispectrum estimator of a single-point time series, - block-averaged with the same window and blocking as BMD -- the quantity - the paper compares the mode bispectrum against in its noise test. - ''' - win = np.hanning(n_dft + 1)[:-1] - win_weight = 1.0 / win.mean() - n_blocks = len(q_x0) // n_dft - q_c = q_x0 - q_x0.mean() - blocks = np.stack([q_c[i * n_dft:(i + 1) * n_dft] for i in range(n_blocks)]) - q_hat = np.fft.fftshift( - np.fft.fft(win[None, :] * blocks, axis=1) * win_weight / n_dft, axes=1) - f_idx = np.rint(np.fft.fftshift(np.fft.fftfreq(n_dft, dt)) * n_dft).astype(int) - row = lambda k: int(np.searchsorted(f_idx, k)) - B = np.full((m, m), np.nan) - for i in range(m): - for j in range(i + 1): - if i + j < n_dft // 2: - B[i, j] = np.abs(np.mean( - q_hat[:, row(i)] * q_hat[:, row(j)] * np.conj(q_hat[:, row(i + j)]))) - return B - - def _bispectrum_grid(bmd, m=40): '''``|lambda_1|`` re-indexed onto a dense ``(k, l)`` grid, NaN elsewhere.''' t = bmd.triads @@ -198,7 +118,7 @@ def test_classical_bispectrum_matches_mode_bispectrum_without_noise(): bmd, x, k = _case(**TRIAD) L_grid = _bispectrum_grid(bmd) q, _, _ = surrogate_waves(TRIAD['freqs'], seed=0) - B = classical_bispectrum(q[:, 0, 0], 128, 1.0, 40) + B = classical_bispectrum(q[:, 0, 0], 128, 40) assert (np.unravel_index(np.nanargmax(L_grid), L_grid.shape) == np.unravel_index(np.nanargmax(B), B.shape)) @@ -213,7 +133,7 @@ def test_triad_survives_unit_snr_noise(): dominate over the rest of the plane -- the paper reports "no significant side peaks". ''' - bmd, x, k = _case(name='noise', freqs=TRIAD['freqs'], snr=1.0) + bmd, x, k = _case(**NOISE) t = bmd.triads vals = np.abs(bmd.L[t.f1_idx, t.f2_idx]) i = int(np.argmax(vals)) @@ -229,7 +149,7 @@ def test_triad_peak_height_is_stable_with_unit_snr_noise(): eigenvalue should remain on the clean-signal scale. ''' clean_bmd, _, _ = _case(**TRIAD) - noisy_bmd, _, _ = _case(name='noise', freqs=TRIAD['freqs'], snr=1.0) + noisy_bmd, _, _ = _case(**NOISE) clean_t = clean_bmd.triads noisy_t = noisy_bmd.triads @@ -263,122 +183,3 @@ def overlap(p, ref): assert overlap(psi_sum, np.exp(-1j * k[2] * x)) == pytest.approx(1.0, abs=1e-3) assert overlap(psi_prod, np.exp(-1j * (k[0] + k[1]) * x)) == pytest.approx(1.0, abs=1e-3) - - -# --------------------------------------------------------------------------- -# reference figures -# --------------------------------------------------------------------------- - -@pytest.mark.slow -def test_reference_figures(): - ''' - Render figures next to Schmidt (2020)'s "hypothesis testing" figures for - visual comparison. Only existence/non-emptiness is asserted here; the - scientific content is covered by the tests above. - ''' - import matplotlib - matplotlib.use('Agg') - import matplotlib.pyplot as plt - - os.makedirs(FIGURES_DIR, exist_ok=True) - - # -- figure 1: 3 rows (nonres / triad / quartet) x 2 columns ------------- - titles = { - 'nonres': r'$f_1 \pm f_2 \pm f_3 \neq 0$ (no triad)', - 'triad': r'$f_1 + f_2 = f_3$ (triad)', - 'quartet': (r'$f_1+f_2+f_3=f_4$, $f_k\pm f_l\pm f_m\neq 0$' - '\n(quartet, no triad)'), - } - colors = ['tab:blue', 'tab:red', 'tab:green', 'tab:purple'] - fig = plt.figure(figsize=(9, 12)) - for row, case in enumerate((NONRES, TRIAD, QUARTET)): - bmd, x, k = _case(**case) - q, _, _ = surrogate_waves(case['freqs'], seed=0) - - ax_a = fig.add_subplot(3, 2, 2 * row + 1) - freq, A = amplitude_spectrum(q[:, 0, 0], 128, 1.0) - pos = freq >= 0 - ax_a.plot(freq[pos], A[pos], 'k') - for j, f in enumerate(case['freqs']): - ax_a.axvline(f, color=colors[j], lw=1) - ax_a.set_xlim(0, 0.5) - ax_a.set_ylim(0, 1.05) - ax_a.set_xlabel('$f$') - ax_a.set_ylabel('$A$') - ax_a.set_title(titles[case['name']], fontsize=9) - - ax_b = fig.add_subplot(3, 2, 2 * row + 2, projection='3d') - t = bmd.triads - vals = np.abs(bmd.L[t.f1_idx, t.f2_idx]) - # plot_trisurf colours by the *data* range, not by set_zlim, so a - # panel that is flat relative to the z-axis (e.g. the non-resonant and - # quartet cases, at <1% of the 0-0.05 range) would otherwise be - # painted with the full colormap and read as structured; pin vmin/vmax - # to the z-limits so colour and height agree, as MATLAB's fixed caxis - # does in the published figure. - ax_b.plot_trisurf(t.f1, t.f2, vals, cmap='viridis', linewidth=0.1, - vmin=0, vmax=0.05) - ax_b.set_zlim(0, 0.05) - ax_b.set_xlabel('$f_1$') - ax_b.set_ylabel('$f_2$') - ax_b.set_zlabel(r'$|\lambda_1|$') - ax_b.set_xticks([0, 0.2, 0.4]) - ax_b.set_yticks([0, 0.1, 0.2]) - fig.tight_layout() - out1 = os.path.join(FIGURES_DIR, 'hypothesis_harmonics_row.png') - fig.savefig(out1, dpi=150) - plt.close(fig) - assert os.path.getsize(out1) > 0 - - # -- figure 2: noise case, 3 panels (bispectra as 3-D surfaces, matching - # the style of figure 1's mode-bispectrum panels) ---------------------- - bmd, x, k = _case(name='noise', freqs=TRIAD['freqs'], snr=1.0) - q, _, _ = surrogate_waves(TRIAD['freqs'], seed=0, snr=1.0) - t = bmd.triads - B = classical_bispectrum(q[:, 0, 0], 128, 1.0, 64) - B_vals = B[t.k, t.l] - - fig2 = plt.figure(figsize=(15, 4.5)) - ax0 = fig2.add_subplot(1, 3, 1) - freq, A = amplitude_spectrum(q[:, 0, 0], 128, 1.0) - pos = freq >= 0 - ax0.plot(freq[pos], A[pos], 'k') - ax0.set_xlim(0, 0.5) - ax0.set_title('(a) amplitude spectrum') - ax0.set_xlabel('$f$') - - ax1 = fig2.add_subplot(1, 3, 2, projection='3d') - ax1.plot_trisurf(t.f1, t.f2, B_vals, cmap='viridis', linewidth=0.1, - vmin=0, vmax=0.25) - ax1.set_title('(b) classical bispectrum') - ax1.set_zlim(0, 0.25) - ax1.set_xlabel('$f_1$') - ax1.set_ylabel('$f_2$') - ax1.set_zlabel('$|B|$') - ax1.set_xticks([0, 0.2, 0.4]) - ax1.set_yticks([0, 0.1, 0.2]) - - ax2 = fig2.add_subplot(1, 3, 3, projection='3d') - # the z-limit must come from the data: the noisy mode bispectrum peaks - # around 0.057 here, and - # a fixed limit taken from an unrelated scale (e.g. panel (b)'s, or a - # guessed round number) either flattens the peak into invisibility or - # leaves the panel mostly empty -- neither matches the published figure, - # which shows one clear dominant peak. - vals = np.abs(bmd.L[t.f1_idx, t.f2_idx]) - z_max = float(np.nanmax(vals)) * 1.05 - ax2.plot_trisurf(t.f1, t.f2, vals, cmap='viridis', linewidth=0.1, - vmin=0, vmax=z_max) - ax2.set_title('(c) mode bispectrum') - ax2.set_zlim(0, z_max) - ax2.set_xlabel('$f_1$') - ax2.set_ylabel('$f_2$') - ax2.set_zlabel(r'$|\lambda_1|$') - ax2.set_xticks([0, 0.2, 0.4]) - ax2.set_yticks([0, 0.1, 0.2]) - - fig2.tight_layout() - out2 = os.path.join(FIGURES_DIR, 'hypothesis_noise.png') - fig2.savefig(out2, dpi=150) - plt.close(fig2) - assert os.path.getsize(out2) > 0 diff --git a/tests/test_io_rejects_complex.py b/tests/test_io_rejects_complex.py new file mode 100644 index 0000000..b1ba7de --- /dev/null +++ b/tests/test_io_rejects_complex.py @@ -0,0 +1,18 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +'''Complex data must be rejected, not silently cast to its real part.''' +import numpy as np +import pytest + +from pybmd.bmd.standard import Standard +from pybmd.utils.io import get_data_array + + +def test_complex_data_is_rejected(tmp_path): + data = np.ones((64, 4, 1)) * 1j + with pytest.raises(TypeError, match='real-valued'): + get_data_array(data, xdim=1, nv=1) + params = dict(n_dft=16, time_step=1.0, n_space_dims=1, n_variables=1, + savedir=str(tmp_path)) + with pytest.raises(TypeError, match='real-valued'): + Standard(params=params).fit(data) diff --git a/tests/test_octave_reference.py b/tests/test_octave_reference.py index c195e7c..16d0f63 100644 --- a/tests/test_octave_reference.py +++ b/tests/test_octave_reference.py @@ -21,7 +21,7 @@ matrices, with a dense angular scan plus local refinement as an independent check -- this is where the deviations documented in ``CLAUDE.md`` live. -See ``docs/octave_cross_validation.md`` for the full measured tables. +See ``tests/octave/octave_cross_validation.md`` for the full measured tables. ''' import os import sys @@ -71,7 +71,6 @@ def _small_bmd(small_data, tmp_path, **overrides): w = utils_weights.uniform((small_data['n1'], small_data['n2']), 1, 1.0) bmd = Standard(params=params, weights=w) bmd._initialize(x) - bmd._mode_shape = (*bmd._xshape, bmd._nv) # normally set by fit() return bmd, x, w @@ -179,7 +178,7 @@ def test_tier_a_qhat_and_b_match(small_data, tmp_path, weight_kind): out = _octave_bmd(bmd, x, w, instrumented=True) q_hat_ref, b_all_ref = out['Q_hat'], out['B_all'] - q_hat = bmd._compute_qhat(block_shape=(bmd._nxv,)) + q_hat = bmd._compute_qhat() for f in t.freq_needed: f = int(f) @@ -255,7 +254,7 @@ def test_tier_c_full_dataset_matches_measured_deviation( (regions=[1,2], max_freq_idx=12) on the full cylinder-wake dataset, and pins the measured counts as a regression: 52/169 triads off by >1%, 29/169 by >10%, always an under-estimate. See - docs/octave_cross_validation.md for the full table this comes from. + tests/octave/octave_cross_validation.md for the full table this comes from. ''' import scipy.io mat_path = oref.require_full_dataset() @@ -345,7 +344,7 @@ def test_tier_b_modes_and_energy_transfer_agree(small_data, tmp_path): out = _octave_bmd(bmd, x, w) L_ref, T_ref, P_ref = out['L'], out['T'], out['P'] - q_hat = bmd._compute_qhat(block_shape=(bmd._nxv,)) + q_hat = bmd._compute_qhat() bmd._triad_loop(q_hat) vals_py = np.abs(bmd.L[t.f1_idx, t.f2_idx]) @@ -378,8 +377,6 @@ def _small_cbmd(small_data, tmp_path, x3, **overrides): params.update(overrides) cb = Cross(params=params) cb._initialize(x3) - cb._mode_shape = (*cb._xshape, cb.n_state) # normally set by fit() - cb._weights_tiled = np.tile(cb._weights, (cb.n_state, 1)) return cb @@ -404,7 +401,7 @@ def test_cbmd_tier_a_matches_reference(small_data, tmp_path, case): cb = _small_cbmd(small_data, tmp_path, x3) out = _octave_cbmd(cb, x3, instrumented=True) b_all_ref = out['B_all'] - q_hat = cb._compute_qhat(block_shape=(cb._nx, cb._nv)) + q_hat = cb._compute_qhat() max_rel = 0.0 for i in range(cb.n_triads):