|
| 1 | +# -*- coding: utf-8 -*- |
| 2 | +""" |
| 3 | +Automated tests for Variational Mode Extraction (VME). |
| 4 | +
|
| 5 | +Covers every public function and method in ``pysdkit._vmd.vme``, including |
| 6 | +helpers ported from MATLAB ``vme.m`` and the File Exchange demos in |
| 7 | +``VME_test_script.m``. |
| 8 | +""" |
| 9 | + |
| 10 | +from __future__ import annotations |
| 11 | + |
| 12 | +import unittest |
| 13 | + |
| 14 | +import numpy as np |
| 15 | + |
| 16 | +from pysdkit import VME, vme |
| 17 | +from pysdkit._vmd.vme import ( |
| 18 | + compactness_kernel, |
| 19 | + crop_mirror, |
| 20 | + ensure_even_length, |
| 21 | + generate_vme_example1, |
| 22 | + generate_vme_example2, |
| 23 | + generate_vme_example3a, |
| 24 | + generate_vme_example3b, |
| 25 | + load_vme_ecg_055m, |
| 26 | + mirror_extend, |
| 27 | + onesided_fft, |
| 28 | + reconstruct_hermitian, |
| 29 | + relative_spectrum_change, |
| 30 | + residual_spectrum, |
| 31 | + spectral_axis, |
| 32 | + spectrum_to_time, |
| 33 | + update_center_frequency, |
| 34 | + update_dual, |
| 35 | + update_mode_spectrum, |
| 36 | +) |
| 37 | + |
| 38 | + |
| 39 | +def _corr(a: np.ndarray, b: np.ndarray) -> float: |
| 40 | + a = np.asarray(a, dtype=float).ravel() |
| 41 | + b = np.asarray(b, dtype=float).ravel() |
| 42 | + a = a - a.mean() |
| 43 | + b = b - b.mean() |
| 44 | + denom = np.linalg.norm(a) * np.linalg.norm(b) |
| 45 | + if denom == 0.0: |
| 46 | + return 0.0 |
| 47 | + return float(np.dot(a, b) / denom) |
| 48 | + |
| 49 | + |
| 50 | +class EnsureEvenLengthTest(unittest.TestCase): |
| 51 | + def test_even_passthrough(self) -> None: |
| 52 | + x = np.arange(8, dtype=float) |
| 53 | + np.testing.assert_array_equal(ensure_even_length(x), x) |
| 54 | + |
| 55 | + def test_odd_drops_last(self) -> None: |
| 56 | + x = np.arange(9, dtype=float) |
| 57 | + np.testing.assert_array_equal(ensure_even_length(x), x[:-1]) |
| 58 | + |
| 59 | + def test_too_short(self) -> None: |
| 60 | + with self.assertRaises(ValueError): |
| 61 | + ensure_even_length(np.array([1.0])) |
| 62 | + |
| 63 | + |
| 64 | +class MirrorExtendTest(unittest.TestCase): |
| 65 | + def test_length_and_halves(self) -> None: |
| 66 | + x = np.arange(10, dtype=float) |
| 67 | + ext = mirror_extend(x) |
| 68 | + self.assertEqual(ext.size, 20) |
| 69 | + np.testing.assert_array_equal(ext[:5], x[:5][::-1]) |
| 70 | + np.testing.assert_array_equal(ext[5:15], x) |
| 71 | + np.testing.assert_array_equal(ext[15:], x[5:][::-1]) |
| 72 | + |
| 73 | + def test_crop_inverts_extend(self) -> None: |
| 74 | + x = np.linspace(-1.0, 1.0, 16) |
| 75 | + np.testing.assert_allclose(crop_mirror(mirror_extend(x)), x) |
| 76 | + |
| 77 | + def test_odd_input_then_crop(self) -> None: |
| 78 | + x = np.arange(11, dtype=float) |
| 79 | + cropped = crop_mirror(mirror_extend(x)) |
| 80 | + np.testing.assert_array_equal(cropped, x[:-1]) |
| 81 | + |
| 82 | + |
| 83 | +class SpectralHelpersTest(unittest.TestCase): |
| 84 | + def test_spectral_axis_matches_matlab(self) -> None: |
| 85 | + n = 8 |
| 86 | + t = np.arange(1, n + 1, dtype=float) / n |
| 87 | + expected = t - 0.5 - 1.0 / n |
| 88 | + np.testing.assert_allclose(spectral_axis(n), expected) |
| 89 | + self.assertAlmostEqual(spectral_axis(n)[n // 2], 0.0, places=12) |
| 90 | + |
| 91 | + def test_onesided_zeros_negative_bins(self) -> None: |
| 92 | + rng = np.random.default_rng(0) |
| 93 | + x = rng.normal(size=32) |
| 94 | + spec = onesided_fft(x) |
| 95 | + self.assertEqual(spec.dtype, np.complex128) |
| 96 | + np.testing.assert_array_equal(spec[: spec.size // 2], 0.0) |
| 97 | + self.assertTrue(np.any(np.abs(spec[spec.size // 2 :]) > 0.0)) |
| 98 | + |
| 99 | + def test_compactness_kernel_at_center(self) -> None: |
| 100 | + omega = np.array([-0.25, 0.0, 0.25]) |
| 101 | + kernel = compactness_kernel(omega, omega_d=0.0, alpha=10.0) |
| 102 | + self.assertEqual(kernel.shape, omega.shape) |
| 103 | + self.assertAlmostEqual(float(kernel[1]), 0.0) |
| 104 | + self.assertTrue(np.all(kernel >= 0.0)) |
| 105 | + |
| 106 | + |
| 107 | +class UpdateStepsTest(unittest.TestCase): |
| 108 | + def setUp(self) -> None: |
| 109 | + self.n = 16 |
| 110 | + self.omega = spectral_axis(self.n) |
| 111 | + self.f_hat = onesided_fft(np.cos(2.0 * np.pi * np.arange(self.n) / 8.0)) |
| 112 | + self.u0 = np.zeros(self.n, dtype=np.complex128) |
| 113 | + self.dual0 = np.zeros(self.n, dtype=np.complex128) |
| 114 | + self.alpha = 50.0 |
| 115 | + |
| 116 | + def test_mode_spectrum_finite_and_onesided(self) -> None: |
| 117 | + u1 = update_mode_spectrum( |
| 118 | + self.f_hat, self.u0, self.dual0, self.omega, 0.1, self.alpha |
| 119 | + ) |
| 120 | + self.assertEqual(u1.dtype, np.complex128) |
| 121 | + self.assertTrue(np.all(np.isfinite(u1))) |
| 122 | + np.testing.assert_array_equal(u1[: self.n // 2], 0.0) |
| 123 | + |
| 124 | + def test_residual_zero_when_u_equals_f(self) -> None: |
| 125 | + residual = residual_spectrum( |
| 126 | + self.f_hat, self.f_hat, self.omega, 0.1, self.alpha |
| 127 | + ) |
| 128 | + np.testing.assert_allclose(residual, 0.0, atol=1e-12) |
| 129 | + |
| 130 | + def test_dual_unchanged_when_tau_zero(self) -> None: |
| 131 | + u1 = update_mode_spectrum( |
| 132 | + self.f_hat, self.u0, self.dual0, self.omega, 0.1, self.alpha |
| 133 | + ) |
| 134 | + dual1 = update_dual( |
| 135 | + self.dual0, self.f_hat, u1, self.omega, 0.1, self.alpha, tau=0.0 |
| 136 | + ) |
| 137 | + np.testing.assert_array_equal(dual1, self.dual0) |
| 138 | + |
| 139 | + def test_dual_steps_when_tau_positive(self) -> None: |
| 140 | + u1 = update_mode_spectrum( |
| 141 | + self.f_hat, self.u0, self.dual0, self.omega, 0.1, self.alpha |
| 142 | + ) |
| 143 | + dual1 = update_dual( |
| 144 | + self.dual0, self.f_hat, u1, self.omega, 0.1, self.alpha, tau=0.5 |
| 145 | + ) |
| 146 | + self.assertFalse(np.allclose(dual1, self.dual0)) |
| 147 | + |
| 148 | + def test_center_frequency_in_unit_interval(self) -> None: |
| 149 | + u1 = update_mode_spectrum( |
| 150 | + self.f_hat, self.u0, self.dual0, self.omega, 0.1, self.alpha |
| 151 | + ) |
| 152 | + omega = update_center_frequency(u1, self.omega, previous=0.1) |
| 153 | + self.assertGreaterEqual(omega, -0.5) |
| 154 | + self.assertLessEqual(omega, 0.5) |
| 155 | + |
| 156 | + def test_center_frequency_zero_power_keeps_previous(self) -> None: |
| 157 | + zeros = np.zeros(self.n, dtype=np.complex128) |
| 158 | + omega = update_center_frequency(zeros, self.omega, previous=0.17) |
| 159 | + self.assertAlmostEqual(omega, 0.17) |
| 160 | + |
| 161 | + def test_relative_change_zero_for_identical(self) -> None: |
| 162 | + x = np.ones(8, dtype=np.complex128) |
| 163 | + self.assertAlmostEqual(relative_spectrum_change(x, x), np.finfo(float).eps) |
| 164 | + |
| 165 | + |
| 166 | +class ReconstructTest(unittest.TestCase): |
| 167 | + def test_hermitian_symmetry_after_overwrite(self) -> None: |
| 168 | + t = 8 |
| 169 | + pos = np.zeros(t, dtype=np.complex128) |
| 170 | + pos[4:] = np.array([1, 2 + 1j, 3 - 2j, 4 + 0.5j]) |
| 171 | + u_hat = reconstruct_hermitian(pos) |
| 172 | + self.assertEqual(u_hat.size, t) |
| 173 | + np.testing.assert_allclose(u_hat[0], np.conj(u_hat[-1])) |
| 174 | + # MATLAB overwrites Nyquist bin with conj(pos[half]) |
| 175 | + np.testing.assert_allclose(u_hat[4], np.conj(pos[4])) |
| 176 | + |
| 177 | + def test_spectrum_to_time_real(self) -> None: |
| 178 | + x = np.cos(2.0 * np.pi * np.arange(32) / 8.0) |
| 179 | + spec = np.fft.fftshift(np.fft.fft(x)) |
| 180 | + rec = spectrum_to_time(spec) |
| 181 | + np.testing.assert_allclose(rec, x, atol=1e-10) |
| 182 | + |
| 183 | + |
| 184 | +class DemoGeneratorsTest(unittest.TestCase): |
| 185 | + def test_example_shapes(self) -> None: |
| 186 | + for factory in ( |
| 187 | + generate_vme_example1, |
| 188 | + generate_vme_example2, |
| 189 | + generate_vme_example3a, |
| 190 | + generate_vme_example3b, |
| 191 | + ): |
| 192 | + demo = factory(n_samples=1000, fs=1000.0) |
| 193 | + self.assertEqual(demo["signal"].shape, (1000,)) |
| 194 | + self.assertEqual(demo["reference"].shape, (1000,)) |
| 195 | + self.assertEqual(demo["t"].shape, (1000,)) |
| 196 | + self.assertEqual(demo["fs"], 1000.0) |
| 197 | + |
| 198 | + def test_example2_three_tones(self) -> None: |
| 199 | + demo = generate_vme_example2() |
| 200 | + self.assertGreater(np.std(demo["signal"]), np.std(demo["reference"])) |
| 201 | + |
| 202 | + |
| 203 | +class LoadEcgTest(unittest.TestCase): |
| 204 | + def test_packaged_mimic_055m(self) -> None: |
| 205 | + demo = load_vme_ecg_055m() |
| 206 | + self.assertEqual(demo["val"].shape, (7, 7500)) |
| 207 | + self.assertEqual(demo["ecg"].shape, (7500,)) |
| 208 | + self.assertEqual(demo["respiration"].shape, (7500,)) |
| 209 | + self.assertEqual(demo["fs"], 125.0) |
| 210 | + self.assertEqual(demo["t"].size, 7500) |
| 211 | + np.testing.assert_array_equal(demo["ecg"], demo["val"][0]) |
| 212 | + np.testing.assert_array_equal(demo["respiration"], demo["val"][-1]) |
| 213 | + self.assertAlmostEqual(demo["t"][-1], (7500 - 1) / 125.0) |
| 214 | + |
| 215 | + |
| 216 | +class VMEClassTest(unittest.TestCase): |
| 217 | + def test_str(self) -> None: |
| 218 | + self.assertIn("VME", str(VME())) |
| 219 | + |
| 220 | + def test_invalid_params(self) -> None: |
| 221 | + with self.assertRaises(ValueError): |
| 222 | + VME(alpha=-1.0) |
| 223 | + with self.assertRaises(ValueError): |
| 224 | + VME(fs=0.0) |
| 225 | + with self.assertRaises(ValueError): |
| 226 | + VME(max_iter=1) |
| 227 | + with self.assertRaises(ValueError): |
| 228 | + VME(tol=0.0) |
| 229 | + |
| 230 | + def test_call_matches_fit_transform(self) -> None: |
| 231 | + demo = generate_vme_example2() |
| 232 | + a = VME(alpha=20000.0, omega_init=10.0, fs=1000.0)(demo["signal"]) |
| 233 | + b = VME(alpha=20000.0, omega_init=10.0, fs=1000.0).fit_transform(demo["signal"]) |
| 234 | + np.testing.assert_allclose(a, b) |
| 235 | + |
| 236 | + def test_return_all_shapes(self) -> None: |
| 237 | + demo = generate_vme_example2() |
| 238 | + extractor = VME(alpha=20000.0, omega_init=10.0, fs=1000.0, max_iter=80) |
| 239 | + u, u_hat, omega = extractor.fit_transform(demo["signal"], return_all=True) |
| 240 | + self.assertEqual(u.shape, demo["signal"].shape) |
| 241 | + self.assertEqual(u_hat.shape, demo["signal"].shape) |
| 242 | + self.assertGreaterEqual(omega.size, 2) |
| 243 | + self.assertEqual(extractor.n_iter, omega.size - 1) |
| 244 | + self.assertAlmostEqual(extractor.omega, float(omega[-1])) |
| 245 | + self.assertTrue(np.all(np.isfinite(u))) |
| 246 | + self.assertTrue(np.all(np.isfinite(omega))) |
| 247 | + |
| 248 | + def test_functional_interface(self) -> None: |
| 249 | + demo = generate_vme_example1() |
| 250 | + u_cls = VME(alpha=20000.0, omega_init=0.0, fs=1000.0).fit_transform( |
| 251 | + demo["signal"] |
| 252 | + ) |
| 253 | + u_fn = vme(demo["signal"], alpha=20000.0, omega_init=0.0, fs=1000.0) |
| 254 | + np.testing.assert_allclose(u_cls, u_fn) |
| 255 | + |
| 256 | + def test_odd_length_truncated(self) -> None: |
| 257 | + demo = generate_vme_example2(n_samples=1001, fs=1000.0) |
| 258 | + u = VME(alpha=5000.0, omega_init=10.0, fs=1000.0, max_iter=40).fit_transform( |
| 259 | + demo["signal"] |
| 260 | + ) |
| 261 | + self.assertEqual(u.size, 1000) |
| 262 | + |
| 263 | + def test_stores_attributes(self) -> None: |
| 264 | + demo = generate_vme_example1() |
| 265 | + extractor = VME(alpha=20000.0, omega_init=0.0, fs=1000.0, max_iter=50) |
| 266 | + u = extractor.fit_transform(demo["signal"]) |
| 267 | + np.testing.assert_array_equal(extractor.u, u) |
| 268 | + self.assertIsNotNone(extractor.u_hat) |
| 269 | + self.assertIsNotNone(extractor.omega_hist) |
| 270 | + |
| 271 | + |
| 272 | +class MatlabDemoAccuracyTest(unittest.TestCase): |
| 273 | + """Reproduce ``VME_test_script.m`` and check correlation with the target mode.""" |
| 274 | + |
| 275 | + def test_example1_lowest_mode(self) -> None: |
| 276 | + demo = generate_vme_example1() |
| 277 | + u = VME( |
| 278 | + alpha=20000.0, omega_init=0.0, fs=1000.0, tau=0.0, tol=1e-7 |
| 279 | + ).fit_transform(demo["signal"]) |
| 280 | + self.assertGreater(_corr(u, demo["reference"]), 0.95) |
| 281 | + |
| 282 | + def test_example2_am_tone(self) -> None: |
| 283 | + demo = generate_vme_example2() |
| 284 | + u = VME( |
| 285 | + alpha=20000.0, omega_init=10.0, fs=1000.0, tau=0.0, tol=1e-7 |
| 286 | + ).fit_transform(demo["signal"]) |
| 287 | + self.assertGreater(_corr(u, demo["reference"]), 0.90) |
| 288 | + # Should not lock onto the 2 Hz trend |
| 289 | + c1 = 2.0 * np.cos(4.0 * np.pi * demo["t"]) |
| 290 | + self.assertLess(_corr(u, c1), 0.4) |
| 291 | + |
| 292 | + def test_example3a_chirp(self) -> None: |
| 293 | + demo = generate_vme_example3a() |
| 294 | + u = VME( |
| 295 | + alpha=20000.0, omega_init=6.0, fs=1000.0, tau=0.0, tol=1e-7 |
| 296 | + ).fit_transform(demo["signal"]) |
| 297 | + self.assertGreater(_corr(u, demo["reference"]), 0.85) |
| 298 | + |
| 299 | + def test_example3b_piecewise_tone(self) -> None: |
| 300 | + demo = generate_vme_example3b() |
| 301 | + u = VME( |
| 302 | + alpha=20000.0, omega_init=26.0, fs=1000.0, tau=0.0, tol=1e-7 |
| 303 | + ).fit_transform(demo["signal"]) |
| 304 | + self.assertGreater(_corr(u, demo["reference"]), 0.70) |
| 305 | + |
| 306 | + def test_ecg_derived_respiration(self) -> None: |
| 307 | + rec = load_vme_ecg_055m() |
| 308 | + # Paper Fig. 4 uses a 32 s (4000-sample) window of record 055 |
| 309 | + n = 4000 |
| 310 | + extractor = VME( |
| 311 | + alpha=20000.0, |
| 312 | + omega_init=0.0, |
| 313 | + fs=float(rec["fs"]), |
| 314 | + tau=0.0, |
| 315 | + tol=1e-7, |
| 316 | + ) |
| 317 | + u, _, omega_hist = extractor.fit_transform(rec["ecg"][:n], return_all=True) |
| 318 | + resp = rec["respiration"][:n] |
| 319 | + self.assertEqual(u.size, n) |
| 320 | + self.assertGreater(_corr(u, resp), 0.3) |
| 321 | + omega_hz = float(omega_hist[-1]) * float(rec["fs"]) |
| 322 | + self.assertGreaterEqual(omega_hz, 0.0) |
| 323 | + self.assertLess(omega_hz, 1.0) |
| 324 | + |
| 325 | + |
| 326 | +if __name__ == "__main__": |
| 327 | + unittest.main() |
0 commit comments