test_peak_finding.py 35 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913
  1. import copy
  2. import numpy as np
  3. import pytest
  4. from pytest import raises, warns
  5. from scipy._lib._array_api import xp_assert_close, xp_assert_equal
  6. from scipy.signal._peak_finding import (
  7. argrelmax,
  8. argrelmin,
  9. peak_prominences,
  10. peak_widths,
  11. _unpack_condition_args,
  12. find_peaks,
  13. find_peaks_cwt,
  14. _identify_ridge_lines
  15. )
  16. from scipy.signal.windows import gaussian
  17. from scipy.signal._peak_finding_utils import _local_maxima_1d, PeakPropertyWarning
  18. def _gen_gaussians(center_locs, sigmas, total_length):
  19. xdata = np.arange(0, total_length).astype(float)
  20. out_data = np.zeros(total_length, dtype=float)
  21. for ind, sigma in enumerate(sigmas):
  22. tmp = (xdata - center_locs[ind]) / sigma
  23. out_data += np.exp(-(tmp**2))
  24. return out_data
  25. def _gen_gaussians_even(sigmas, total_length):
  26. num_peaks = len(sigmas)
  27. delta = total_length / (num_peaks + 1)
  28. center_locs = np.linspace(delta, total_length - delta, num=num_peaks).astype(int)
  29. out_data = _gen_gaussians(center_locs, sigmas, total_length)
  30. return out_data, center_locs
  31. def _gen_ridge_line(start_locs, max_locs, length, distances, gaps):
  32. """
  33. Generate coordinates for a ridge line.
  34. Will be a series of coordinates, starting a start_loc (length 2).
  35. The maximum distance between any adjacent columns will be
  36. `max_distance`, the max distance between adjacent rows
  37. will be `map_gap'.
  38. `max_locs` should be the size of the intended matrix. The
  39. ending coordinates are guaranteed to be less than `max_locs`,
  40. although they may not approach `max_locs` at all.
  41. """
  42. def keep_bounds(num, max_val):
  43. out = max(num, 0)
  44. out = min(out, max_val)
  45. return out
  46. gaps = copy.deepcopy(gaps)
  47. distances = copy.deepcopy(distances)
  48. locs = np.zeros([length, 2], dtype=int)
  49. locs[0, :] = start_locs
  50. total_length = max_locs[0] - start_locs[0] - sum(gaps)
  51. if total_length < length:
  52. raise ValueError('Cannot generate ridge line according to constraints')
  53. dist_int = length / len(distances) - 1
  54. gap_int = length / len(gaps) - 1
  55. for ind in range(1, length):
  56. nextcol = locs[ind - 1, 1]
  57. nextrow = locs[ind - 1, 0] + 1
  58. if (ind % dist_int == 0) and (len(distances) > 0):
  59. nextcol += ((-1)**ind)*distances.pop()
  60. if (ind % gap_int == 0) and (len(gaps) > 0):
  61. nextrow += gaps.pop()
  62. nextrow = keep_bounds(nextrow, max_locs[0])
  63. nextcol = keep_bounds(nextcol, max_locs[1])
  64. locs[ind, :] = [nextrow, nextcol]
  65. return [locs[:, 0], locs[:, 1]]
  66. class TestLocalMaxima1d:
  67. def test_empty(self):
  68. """Test with empty signal."""
  69. x = np.array([], dtype=np.float64)
  70. for array in _local_maxima_1d(x):
  71. xp_assert_equal(array, np.array([]), check_dtype=False)
  72. assert array.base is None
  73. def test_linear(self):
  74. """Test with linear signal."""
  75. x = np.linspace(0, 100)
  76. for array in _local_maxima_1d(x):
  77. xp_assert_equal(array, np.array([], dtype=np.intp))
  78. assert array.base is None
  79. def test_simple(self):
  80. """Test with simple signal."""
  81. x = np.linspace(-10, 10, 50)
  82. x[2::3] += 1
  83. expected = np.arange(2, 50, 3, dtype=np.intp)
  84. for array in _local_maxima_1d(x):
  85. # For plateaus of size 1, the edges are identical with the
  86. # midpoints
  87. xp_assert_equal(array, expected, check_dtype=False)
  88. assert array.base is None
  89. def test_flat_maxima(self):
  90. """Test if flat maxima are detected correctly."""
  91. x = np.array([-1.3, 0, 1, 0, 2, 2, 0, 3, 3, 3, 2.99, 4, 4, 4, 4, -10,
  92. -5, -5, -5, -5, -5, -10])
  93. midpoints, left_edges, right_edges = _local_maxima_1d(x)
  94. xp_assert_equal(midpoints, np.array([2, 4, 8, 12, 18]), check_dtype=False)
  95. xp_assert_equal(left_edges, np.array([2, 4, 7, 11, 16]), check_dtype=False)
  96. xp_assert_equal(right_edges, np.array([2, 5, 9, 14, 20]), check_dtype=False)
  97. @pytest.mark.parametrize('x', [
  98. np.array([1., 0, 2]),
  99. np.array([3., 3, 0, 4, 4]),
  100. np.array([5., 5, 5, 0, 6, 6, 6]),
  101. ])
  102. def test_signal_edges(self, x):
  103. """Test if behavior on signal edges is correct."""
  104. for array in _local_maxima_1d(x):
  105. xp_assert_equal(array, np.array([], dtype=np.intp))
  106. assert array.base is None
  107. def test_exceptions(self):
  108. """Test input validation and raised exceptions."""
  109. with raises(ValueError, match="wrong number of dimensions"):
  110. _local_maxima_1d(np.ones((1, 1)))
  111. with raises(ValueError, match="expected 'const float64_t'"):
  112. _local_maxima_1d(np.ones(1, dtype=int))
  113. with raises(TypeError, match="list"):
  114. _local_maxima_1d([1., 2.])
  115. with raises(TypeError, match="'x' must not be None"):
  116. _local_maxima_1d(None)
  117. class TestRidgeLines:
  118. def test_empty(self):
  119. test_matr = np.zeros([20, 100])
  120. lines = _identify_ridge_lines(test_matr, np.full(20, 2), 1)
  121. assert len(lines) == 0
  122. def test_minimal(self):
  123. test_matr = np.zeros([20, 100])
  124. test_matr[0, 10] = 1
  125. lines = _identify_ridge_lines(test_matr, np.full(20, 2), 1)
  126. assert len(lines) == 1
  127. test_matr = np.zeros([20, 100])
  128. test_matr[0:2, 10] = 1
  129. lines = _identify_ridge_lines(test_matr, np.full(20, 2), 1)
  130. assert len(lines) == 1
  131. def test_single_pass(self):
  132. distances = [0, 1, 2, 5]
  133. gaps = [0, 1, 2, 0, 1]
  134. test_matr = np.zeros([20, 50]) + 1e-12
  135. length = 12
  136. line = _gen_ridge_line([0, 25], test_matr.shape, length, distances, gaps)
  137. test_matr[line[0], line[1]] = 1
  138. max_distances = np.full(20, max(distances))
  139. identified_lines = _identify_ridge_lines(test_matr,
  140. max_distances,
  141. max(gaps) + 1)
  142. assert len(identified_lines) == 1
  143. for iline_, line_ in zip(identified_lines[0], line):
  144. xp_assert_equal(iline_, line_, check_dtype=False)
  145. def test_single_bigdist(self):
  146. distances = [0, 1, 2, 5]
  147. gaps = [0, 1, 2, 4]
  148. test_matr = np.zeros([20, 50])
  149. length = 12
  150. line = _gen_ridge_line([0, 25], test_matr.shape, length, distances, gaps)
  151. test_matr[line[0], line[1]] = 1
  152. max_dist = 3
  153. max_distances = np.full(20, max_dist)
  154. #This should get 2 lines, since the distance is too large
  155. identified_lines = _identify_ridge_lines(test_matr,
  156. max_distances,
  157. max(gaps) + 1)
  158. assert len(identified_lines) == 2
  159. for iline in identified_lines:
  160. adists = np.diff(iline[1])
  161. np.testing.assert_array_less(np.abs(adists), max_dist)
  162. agaps = np.diff(iline[0])
  163. np.testing.assert_array_less(np.abs(agaps), max(gaps) + 0.1)
  164. def test_single_biggap(self):
  165. distances = [0, 1, 2, 5]
  166. max_gap = 3
  167. gaps = [0, 4, 2, 1]
  168. test_matr = np.zeros([20, 50])
  169. length = 12
  170. line = _gen_ridge_line([0, 25], test_matr.shape, length, distances, gaps)
  171. test_matr[line[0], line[1]] = 1
  172. max_dist = 6
  173. max_distances = np.full(20, max_dist)
  174. #This should get 2 lines, since the gap is too large
  175. identified_lines = _identify_ridge_lines(test_matr, max_distances, max_gap)
  176. assert len(identified_lines) == 2
  177. for iline in identified_lines:
  178. adists = np.diff(iline[1])
  179. np.testing.assert_array_less(np.abs(adists), max_dist)
  180. agaps = np.diff(iline[0])
  181. np.testing.assert_array_less(np.abs(agaps), max(gaps) + 0.1)
  182. def test_single_biggaps(self):
  183. distances = [0]
  184. max_gap = 1
  185. gaps = [3, 6]
  186. test_matr = np.zeros([50, 50])
  187. length = 30
  188. line = _gen_ridge_line([0, 25], test_matr.shape, length, distances, gaps)
  189. test_matr[line[0], line[1]] = 1
  190. max_dist = 1
  191. max_distances = np.full(50, max_dist)
  192. #This should get 3 lines, since the gaps are too large
  193. identified_lines = _identify_ridge_lines(test_matr, max_distances, max_gap)
  194. assert len(identified_lines) == 3
  195. for iline in identified_lines:
  196. adists = np.diff(iline[1])
  197. np.testing.assert_array_less(np.abs(adists), max_dist)
  198. agaps = np.diff(iline[0])
  199. np.testing.assert_array_less(np.abs(agaps), max(gaps) + 0.1)
  200. class TestArgrel:
  201. def test_empty(self):
  202. # Regression test for gh-2832.
  203. # When there are no relative extrema, make sure that
  204. # the number of empty arrays returned matches the
  205. # dimension of the input.
  206. empty_array = np.array([], dtype=int)
  207. z1 = np.zeros(5)
  208. i = argrelmin(z1)
  209. xp_assert_equal(len(i), 1)
  210. xp_assert_equal(i[0], empty_array, check_dtype=False)
  211. z2 = np.zeros((3, 5))
  212. row, col = argrelmin(z2, axis=0)
  213. xp_assert_equal(row, empty_array, check_dtype=False)
  214. xp_assert_equal(col, empty_array, check_dtype=False)
  215. row, col = argrelmin(z2, axis=1)
  216. xp_assert_equal(row, empty_array, check_dtype=False)
  217. xp_assert_equal(col, empty_array, check_dtype=False)
  218. def test_basic(self):
  219. # Note: the docstrings for the argrel{min,max,extrema} functions
  220. # do not give a guarantee of the order of the indices, so we'll
  221. # sort them before testing.
  222. x = np.array([[1, 2, 2, 3, 2],
  223. [2, 1, 2, 2, 3],
  224. [3, 2, 1, 2, 2],
  225. [2, 3, 2, 1, 2],
  226. [1, 2, 3, 2, 1]])
  227. row, col = argrelmax(x, axis=0)
  228. order = np.argsort(row)
  229. xp_assert_equal(row[order], [1, 2, 3], check_dtype=False)
  230. xp_assert_equal(col[order], [4, 0, 1], check_dtype=False)
  231. row, col = argrelmax(x, axis=1)
  232. order = np.argsort(row)
  233. xp_assert_equal(row[order], [0, 3, 4], check_dtype=False)
  234. xp_assert_equal(col[order], [3, 1, 2], check_dtype=False)
  235. row, col = argrelmin(x, axis=0)
  236. order = np.argsort(row)
  237. xp_assert_equal(row[order], [1, 2, 3], check_dtype=False)
  238. xp_assert_equal(col[order], [1, 2, 3], check_dtype=False)
  239. row, col = argrelmin(x, axis=1)
  240. order = np.argsort(row)
  241. xp_assert_equal(row[order], [1, 2, 3], check_dtype=False)
  242. xp_assert_equal(col[order], [1, 2, 3], check_dtype=False)
  243. def test_highorder(self):
  244. order = 2
  245. sigmas = [1.0, 2.0, 10.0, 5.0, 15.0]
  246. test_data, act_locs = _gen_gaussians_even(sigmas, 500)
  247. test_data[act_locs + order] = test_data[act_locs]*0.99999
  248. test_data[act_locs - order] = test_data[act_locs]*0.99999
  249. rel_max_locs = argrelmax(test_data, order=order, mode='clip')[0]
  250. assert len(rel_max_locs) == len(act_locs)
  251. assert (rel_max_locs == act_locs).all()
  252. def test_2d_gaussians(self):
  253. sigmas = [1.0, 2.0, 10.0]
  254. test_data, act_locs = _gen_gaussians_even(sigmas, 100)
  255. rot_factor = 20
  256. rot_range = np.arange(0, len(test_data)) - rot_factor
  257. test_data_2 = np.vstack([test_data, test_data[rot_range]])
  258. rel_max_rows, rel_max_cols = argrelmax(test_data_2, axis=1, order=1)
  259. for rw in range(0, test_data_2.shape[0]):
  260. inds = (rel_max_rows == rw)
  261. assert len(rel_max_cols[inds]) == len(act_locs)
  262. assert (act_locs == (rel_max_cols[inds] - rot_factor*rw)).all()
  263. class TestPeakProminences:
  264. def test_empty(self):
  265. """
  266. Test if an empty array is returned if no peaks are provided.
  267. """
  268. out = peak_prominences([1, 2, 3], [])
  269. for arr, dtype in zip(out, [np.float64, np.intp, np.intp]):
  270. assert arr.size == 0
  271. assert arr.dtype == dtype
  272. out = peak_prominences([], [])
  273. for arr, dtype in zip(out, [np.float64, np.intp, np.intp]):
  274. assert arr.size == 0
  275. assert arr.dtype == dtype
  276. def test_basic(self):
  277. """
  278. Test if height of prominences is correctly calculated in signal with
  279. rising baseline (peak widths are 1 sample).
  280. """
  281. # Prepare basic signal
  282. x = np.array([-1, 1.2, 1.2, 1, 3.2, 1.3, 2.88, 2.1])
  283. peaks = np.array([1, 2, 4, 6])
  284. lbases = np.array([0, 0, 0, 5])
  285. rbases = np.array([3, 3, 5, 7])
  286. proms = x[peaks] - np.max([x[lbases], x[rbases]], axis=0)
  287. # Test if calculation matches handcrafted result
  288. out = peak_prominences(x, peaks)
  289. xp_assert_equal(out[0], proms, check_dtype=False)
  290. xp_assert_equal(out[1], lbases, check_dtype=False)
  291. xp_assert_equal(out[2], rbases, check_dtype=False)
  292. def test_edge_cases(self):
  293. """
  294. Test edge cases.
  295. """
  296. # Peaks have same height, prominence and bases
  297. x = [0, 2, 1, 2, 1, 2, 0]
  298. peaks = [1, 3, 5]
  299. proms, lbases, rbases = peak_prominences(x, peaks)
  300. xp_assert_equal(proms, np.asarray([2.0, 2, 2]), check_dtype=False)
  301. xp_assert_equal(lbases, [0, 0, 0], check_dtype=False)
  302. xp_assert_equal(rbases, [6, 6, 6], check_dtype=False)
  303. # Peaks have same height & prominence but different bases
  304. x = [0, 1, 0, 1, 0, 1, 0]
  305. peaks = np.array([1, 3, 5])
  306. proms, lbases, rbases = peak_prominences(x, peaks)
  307. xp_assert_equal(proms, np.asarray([1.0, 1, 1]))
  308. xp_assert_equal(lbases, peaks - 1, check_dtype=False)
  309. xp_assert_equal(rbases, peaks + 1, check_dtype=False)
  310. def test_non_contiguous(self):
  311. """
  312. Test with non-C-contiguous input arrays.
  313. """
  314. x = np.repeat([-9, 9, 9, 0, 3, 1], 2)
  315. peaks = np.repeat([1, 2, 4], 2)
  316. proms, lbases, rbases = peak_prominences(x[::2], peaks[::2])
  317. xp_assert_equal(proms, np.asarray([9.0, 9, 2]))
  318. xp_assert_equal(lbases, [0, 0, 3], check_dtype=False)
  319. xp_assert_equal(rbases, [3, 3, 5], check_dtype=False)
  320. def test_wlen(self):
  321. """
  322. Test if wlen actually shrinks the evaluation range correctly.
  323. """
  324. x = [0, 1, 2, 3, 1, 0, -1]
  325. peak = [3]
  326. # Test rounding behavior of wlen
  327. proms = peak_prominences(x, peak)
  328. for prom, val in zip(proms, [3.0, 0, 6]):
  329. assert prom == val
  330. for wlen, i in [(8, 0), (7, 0), (6, 0), (5, 1), (3.2, 1), (3, 2), (1.1, 2)]:
  331. proms = peak_prominences(x, peak, wlen)
  332. for prom, val in zip(proms, [3. - i, 0 + i, 6 - i]):
  333. assert prom == val
  334. def test_exceptions(self):
  335. """
  336. Verify that exceptions and warnings are raised.
  337. """
  338. # x with dimension > 1
  339. with raises(ValueError, match='1-D array'):
  340. peak_prominences([[0, 1, 1, 0]], [1, 2])
  341. # peaks with dimension > 1
  342. with raises(ValueError, match='1-D array'):
  343. peak_prominences([0, 1, 1, 0], [[1, 2]])
  344. # x with dimension < 1
  345. with raises(ValueError, match='1-D array'):
  346. peak_prominences(3, [0,])
  347. # empty x with supplied
  348. with raises(ValueError, match='not a valid index'):
  349. peak_prominences([], [0])
  350. # invalid indices with non-empty x
  351. for p in [-100, -1, 3, 1000]:
  352. with raises(ValueError, match='not a valid index'):
  353. peak_prominences([1, 0, 2], [p])
  354. # peaks is not cast-able to np.intp
  355. with raises(TypeError, match='cannot safely cast'):
  356. peak_prominences([0, 1, 1, 0], [1.1, 2.3])
  357. # wlen < 3
  358. with raises(ValueError, match='wlen'):
  359. peak_prominences(np.arange(10), [3, 5], wlen=1)
  360. def test_warnings(self):
  361. """
  362. Verify that appropriate warnings are raised.
  363. """
  364. msg = "some peaks have a prominence of 0"
  365. for p in [0, 1, 2]:
  366. with warns(PeakPropertyWarning, match=msg):
  367. peak_prominences([1, 0, 2], [p,])
  368. with warns(PeakPropertyWarning, match=msg):
  369. peak_prominences([0, 1, 1, 1, 0], [2], wlen=2)
  370. class TestPeakWidths:
  371. def test_empty(self):
  372. """
  373. Test if an empty array is returned if no peaks are provided.
  374. """
  375. widths = peak_widths([], [])[0]
  376. assert isinstance(widths, np.ndarray)
  377. assert widths.size == 0
  378. widths = peak_widths([1, 2, 3], [])[0]
  379. assert isinstance(widths, np.ndarray)
  380. assert widths.size == 0
  381. out = peak_widths([], [])
  382. for arr in out:
  383. assert isinstance(arr, np.ndarray)
  384. assert arr.size == 0
  385. @pytest.mark.filterwarnings("ignore:some peaks have a width of 0")
  386. def test_basic(self):
  387. """
  388. Test a simple use case with easy to verify results at different relative
  389. heights.
  390. """
  391. x = np.array([1, 0, 1, 2, 1, 0, -1])
  392. prominence = 2
  393. for rel_height, width_true, lip_true, rip_true in [
  394. (0., 0., 3., 3.), # raises warning
  395. (0.25, 1., 2.5, 3.5),
  396. (0.5, 2., 2., 4.),
  397. (0.75, 3., 1.5, 4.5),
  398. (1., 4., 1., 5.),
  399. (2., 5., 1., 6.),
  400. (3., 5., 1., 6.)
  401. ]:
  402. width_calc, height, lip_calc, rip_calc = peak_widths(
  403. x, [3], rel_height)
  404. xp_assert_close(width_calc, np.asarray([width_true]))
  405. xp_assert_close(height, np.asarray([2 - rel_height * prominence]))
  406. xp_assert_close(lip_calc, np.asarray([lip_true]))
  407. xp_assert_close(rip_calc, np.asarray([rip_true]))
  408. def test_non_contiguous(self):
  409. """
  410. Test with non-C-contiguous input arrays.
  411. """
  412. x = np.repeat([0, 100, 50], 4)
  413. peaks = np.repeat([1], 3)
  414. result = peak_widths(x[::4], peaks[::3])
  415. xp_assert_equal(result,
  416. np.asarray([[0.75], [75], [0.75], [1.5]])
  417. )
  418. def test_exceptions(self):
  419. """
  420. Verify that argument validation works as intended.
  421. """
  422. with raises(ValueError, match='1-D array'):
  423. # x with dimension > 1
  424. peak_widths(np.zeros((3, 4)), np.ones(3))
  425. with raises(ValueError, match='1-D array'):
  426. # x with dimension < 1
  427. peak_widths(3, [0])
  428. with raises(ValueError, match='1-D array'):
  429. # peaks with dimension > 1
  430. peak_widths(np.arange(10), np.ones((3, 2), dtype=np.intp))
  431. with raises(ValueError, match='1-D array'):
  432. # peaks with dimension < 1
  433. peak_widths(np.arange(10), 3)
  434. with raises(ValueError, match='not a valid index'):
  435. # peak pos exceeds x.size
  436. peak_widths(np.arange(10), [8, 11])
  437. with raises(ValueError, match='not a valid index'):
  438. # empty x with peaks supplied
  439. peak_widths([], [1, 2])
  440. with raises(TypeError, match='cannot safely cast'):
  441. # peak cannot be safely cast to intp
  442. peak_widths(np.arange(10), [1.1, 2.3])
  443. with raises(ValueError, match='rel_height'):
  444. # rel_height is < 0
  445. peak_widths([0, 1, 0, 1, 0], [1, 3], rel_height=-1)
  446. with raises(TypeError, match='None'):
  447. # prominence data contains None
  448. peak_widths([1, 2, 1], [1], prominence_data=(None, None, None))
  449. def test_warnings(self):
  450. """
  451. Verify that appropriate warnings are raised.
  452. """
  453. msg = "some peaks have a width of 0"
  454. with warns(PeakPropertyWarning, match=msg):
  455. # Case: rel_height is 0
  456. peak_widths([0, 1, 0], [1], rel_height=0)
  457. with warns(PeakPropertyWarning, match=msg):
  458. # Case: prominence is 0 and bases are identical
  459. peak_widths(
  460. [0, 1, 1, 1, 0], [2],
  461. prominence_data=(np.array([0.], np.float64),
  462. np.array([2], np.intp),
  463. np.array([2], np.intp))
  464. )
  465. def test_mismatching_prominence_data(self):
  466. """Test with mismatching peak and / or prominence data."""
  467. x = [0, 1, 0]
  468. peak = [1]
  469. for i, (prominences, left_bases, right_bases) in enumerate([
  470. ((1.,), (-1,), (2,)), # left base not in x
  471. ((1.,), (0,), (3,)), # right base not in x
  472. ((1.,), (2,), (0,)), # swapped bases same as peak
  473. ((1., 1.), (0, 0), (2, 2)), # array shapes don't match peaks
  474. ((1., 1.), (0,), (2,)), # arrays with different shapes
  475. ((1.,), (0, 0), (2,)), # arrays with different shapes
  476. ((1.,), (0,), (2, 2)) # arrays with different shapes
  477. ]):
  478. # Make sure input is matches output of signal.peak_prominences
  479. prominence_data = (np.array(prominences, dtype=np.float64),
  480. np.array(left_bases, dtype=np.intp),
  481. np.array(right_bases, dtype=np.intp))
  482. # Test for correct exception
  483. if i < 3:
  484. match = "prominence data is invalid for peak"
  485. else:
  486. match = "arrays in `prominence_data` must have the same shape"
  487. with raises(ValueError, match=match):
  488. peak_widths(x, peak, prominence_data=prominence_data)
  489. @pytest.mark.filterwarnings("ignore:some peaks have a width of 0")
  490. def test_intersection_rules(self):
  491. """Test if x == eval_height counts as an intersection."""
  492. # Flatt peak with two possible intersection points if evaluated at 1
  493. x = [0, 1, 2, 1, 3, 3, 3, 1, 2, 1, 0]
  494. # relative height is 0 -> width is 0 as well, raises warning
  495. xp_assert_close(peak_widths(x, peaks=[5], rel_height=0),
  496. [(0.,), (3.,), (5.,), (5.,)])
  497. # width_height == x counts as intersection -> nearest 1 is chosen
  498. xp_assert_close(peak_widths(x, peaks=[5], rel_height=2/3),
  499. [(4.,), (1.,), (3.,), (7.,)])
  500. def test_unpack_condition_args():
  501. """
  502. Verify parsing of condition arguments for `scipy.signal.find_peaks` function.
  503. """
  504. x = np.arange(10)
  505. amin_true = x
  506. amax_true = amin_true + 10
  507. peaks = amin_true[1::2]
  508. # Test unpacking with None or interval
  509. assert (None, None) == _unpack_condition_args((None, None), x, peaks)
  510. assert (1, None) == _unpack_condition_args(1, x, peaks)
  511. assert (1, None) == _unpack_condition_args((1, None), x, peaks)
  512. assert (None, 2) == _unpack_condition_args((None, 2), x, peaks)
  513. assert (3., 4.5) == _unpack_condition_args((3., 4.5), x, peaks)
  514. # Test if borders are correctly reduced with `peaks`
  515. amin_calc, amax_calc = _unpack_condition_args((amin_true, amax_true), x, peaks)
  516. xp_assert_equal(amin_calc, amin_true[peaks])
  517. xp_assert_equal(amax_calc, amax_true[peaks])
  518. # Test raises if array borders don't match x
  519. with raises(ValueError, match="array size of lower"):
  520. _unpack_condition_args(amin_true, np.arange(11), peaks)
  521. with raises(ValueError, match="array size of upper"):
  522. _unpack_condition_args((None, amin_true), np.arange(11), peaks)
  523. class TestFindPeaks:
  524. # Keys of optionally returned properties
  525. property_keys = {'peak_heights', 'left_thresholds', 'right_thresholds',
  526. 'prominences', 'left_bases', 'right_bases', 'widths',
  527. 'width_heights', 'left_ips', 'right_ips'}
  528. def test_constant(self):
  529. """
  530. Test behavior for signal without local maxima.
  531. """
  532. open_interval = (None, None)
  533. peaks, props = find_peaks(np.ones(10),
  534. height=open_interval, threshold=open_interval,
  535. prominence=open_interval, width=open_interval)
  536. assert peaks.size == 0
  537. for key in self.property_keys:
  538. assert props[key].size == 0
  539. def test_plateau_size(self):
  540. """
  541. Test plateau size condition for peaks.
  542. """
  543. # Prepare signal with peaks with peak_height == plateau_size
  544. plateau_sizes = np.array([1, 2, 3, 4, 8, 20, 111])
  545. x = np.zeros(plateau_sizes.size * 2 + 1)
  546. x[1::2] = plateau_sizes
  547. repeats = np.ones(x.size, dtype=int)
  548. repeats[1::2] = x[1::2]
  549. x = np.repeat(x, repeats)
  550. # Test full output
  551. peaks, props = find_peaks(x, plateau_size=(None, None))
  552. xp_assert_equal(peaks, [1, 3, 7, 11, 18, 33, 100], check_dtype=False)
  553. xp_assert_equal(props["plateau_sizes"], plateau_sizes, check_dtype=False)
  554. xp_assert_equal(props["left_edges"], peaks - (plateau_sizes - 1) // 2,
  555. check_dtype=False)
  556. xp_assert_equal(props["right_edges"], peaks + plateau_sizes // 2,
  557. check_dtype=False)
  558. # Test conditions
  559. xp_assert_equal(find_peaks(x, plateau_size=4)[0], [11, 18, 33, 100],
  560. check_dtype=False)
  561. xp_assert_equal(find_peaks(x, plateau_size=(None, 3.5))[0], [1, 3, 7],
  562. check_dtype=False)
  563. xp_assert_equal(find_peaks(x, plateau_size=(5, 50))[0], [18, 33],
  564. check_dtype=False)
  565. def test_height_condition(self):
  566. """
  567. Test height condition for peaks.
  568. """
  569. x = (0., 1/3, 0., 2.5, 0, 4., 0)
  570. peaks, props = find_peaks(x, height=(None, None))
  571. xp_assert_equal(peaks, np.array([1, 3, 5]), check_dtype=False)
  572. xp_assert_equal(props['peak_heights'], np.array([1/3, 2.5, 4.]),
  573. check_dtype=False)
  574. xp_assert_equal(find_peaks(x, height=0.5)[0], np.array([3, 5]),
  575. check_dtype=False)
  576. xp_assert_equal(find_peaks(x, height=(None, 3))[0], np.array([1, 3]),
  577. check_dtype=False)
  578. xp_assert_equal(find_peaks(x, height=(2, 3))[0], np.array([3]),
  579. check_dtype=False)
  580. def test_threshold_condition(self):
  581. """
  582. Test threshold condition for peaks.
  583. """
  584. x = (0, 2, 1, 4, -1)
  585. peaks, props = find_peaks(x, threshold=(None, None))
  586. xp_assert_equal(peaks, np.array([1, 3]), check_dtype=False)
  587. xp_assert_equal(props['left_thresholds'], np.array([2.0, 3.0]))
  588. xp_assert_equal(props['right_thresholds'], np.array([1.0, 5.0]))
  589. xp_assert_equal(find_peaks(x, threshold=2)[0], np.array([3]),
  590. check_dtype=False)
  591. xp_assert_equal(find_peaks(x, threshold=3.5)[0], np.array([], dtype=int),
  592. check_dtype=False)
  593. xp_assert_equal(find_peaks(x, threshold=(None, 5))[0], np.array([1, 3]),
  594. check_dtype=False)
  595. xp_assert_equal(find_peaks(x, threshold=(None, 4))[0], np.array([1]),
  596. check_dtype=False)
  597. xp_assert_equal(find_peaks(x, threshold=(2, 4))[0], np.array([], dtype=int),
  598. check_dtype=False)
  599. def test_distance_condition(self):
  600. """
  601. Test distance condition for peaks.
  602. """
  603. # Peaks of different height with constant distance 3
  604. peaks_all = np.arange(1, 21, 3)
  605. x = np.zeros(21)
  606. x[peaks_all] += np.linspace(1, 2, peaks_all.size)
  607. # Test if peaks with "minimal" distance are still selected (distance = 3)
  608. xp_assert_equal(find_peaks(x, distance=3)[0], peaks_all, check_dtype=False)
  609. # Select every second peak (distance > 3)
  610. peaks_subset = find_peaks(x, distance=3.0001)[0]
  611. # Test if peaks_subset is subset of peaks_all
  612. assert np.setdiff1d(peaks_subset, peaks_all, assume_unique=True).size == 0
  613. # Test if every second peak was removed
  614. dfs = np.diff(peaks_subset)
  615. xp_assert_equal(dfs, 6*np.ones_like(dfs))
  616. # Test priority of peak removal
  617. x = [-2, 1, -1, 0, -3]
  618. peaks_subset = find_peaks(x, distance=10)[0] # use distance > x size
  619. assert peaks_subset.size == 1 and peaks_subset[0] == 1
  620. def test_prominence_condition(self):
  621. """
  622. Test prominence condition for peaks.
  623. """
  624. x = np.linspace(0, 10, 100)
  625. peaks_true = np.arange(1, 99, 2)
  626. offset = np.linspace(1, 10, peaks_true.size)
  627. x[peaks_true] += offset
  628. prominences = x[peaks_true] - x[peaks_true + 1]
  629. interval = (3, 9)
  630. keep = np.nonzero(
  631. (interval[0] <= prominences) & (prominences <= interval[1]))
  632. peaks_calc, properties = find_peaks(x, prominence=interval)
  633. xp_assert_equal(peaks_calc, peaks_true[keep], check_dtype=False)
  634. xp_assert_equal(properties['prominences'], prominences[keep], check_dtype=False)
  635. xp_assert_equal(properties['left_bases'],
  636. np.zeros_like(properties['left_bases']))
  637. xp_assert_equal(properties['right_bases'], peaks_true[keep] + 1,
  638. check_dtype=False)
  639. def test_width_condition(self):
  640. """
  641. Test width condition for peaks.
  642. """
  643. x = np.array([1, 0, 1, 2, 1, 0, -1, 4, 0])
  644. peaks, props = find_peaks(x, width=(None, 2), rel_height=0.75)
  645. assert peaks.size == 1
  646. xp_assert_equal(peaks, 7*np.ones_like(peaks))
  647. xp_assert_close(props['widths'], np.asarray([1.35]))
  648. xp_assert_close(props['width_heights'], np.asarray([1.]))
  649. xp_assert_close(props['left_ips'], np.asarray([6.4]))
  650. xp_assert_close(props['right_ips'], np.asarray([7.75]))
  651. def test_properties(self):
  652. """
  653. Test returned properties.
  654. """
  655. open_interval = (None, None)
  656. x = [0, 1, 0, 2, 1.5, 0, 3, 0, 5, 9]
  657. peaks, props = find_peaks(x,
  658. height=open_interval, threshold=open_interval,
  659. prominence=open_interval, width=open_interval)
  660. assert len(props) == len(self.property_keys)
  661. for key in self.property_keys:
  662. assert peaks.size == props[key].size
  663. def test_raises(self):
  664. """
  665. Test exceptions raised by function.
  666. """
  667. with raises(ValueError, match="1-D array"):
  668. find_peaks(np.array(1))
  669. with raises(ValueError, match="1-D array"):
  670. find_peaks(np.ones((2, 2)))
  671. with raises(ValueError, match="distance"):
  672. find_peaks(np.arange(10), distance=-1)
  673. @pytest.mark.filterwarnings("ignore:some peaks have a prominence of 0",
  674. "ignore:some peaks have a width of 0")
  675. def test_wlen_smaller_plateau(self):
  676. """
  677. Test behavior of prominence and width calculation if the given window
  678. length is smaller than a peak's plateau size.
  679. Regression test for gh-9110.
  680. """
  681. peaks, props = find_peaks([0, 1, 1, 1, 0], prominence=(None, None),
  682. width=(None, None), wlen=2)
  683. xp_assert_equal(peaks, 2 * np.ones_like(peaks))
  684. xp_assert_equal(props["prominences"], np.zeros_like(props["prominences"]))
  685. xp_assert_equal(props["widths"], np.zeros_like(props["widths"]))
  686. xp_assert_equal(props["width_heights"], np.ones_like(props["width_heights"]))
  687. for key in ("left_bases", "right_bases", "left_ips", "right_ips"):
  688. xp_assert_equal(props[key], peaks, check_dtype=False)
  689. @pytest.mark.parametrize("kwargs", [
  690. {},
  691. {"distance": 3.0},
  692. {"prominence": (None, None)},
  693. {"width": (None, 2)},
  694. ])
  695. def test_readonly_array(self, kwargs):
  696. """
  697. Test readonly arrays are accepted.
  698. """
  699. x = np.linspace(0, 10, 15)
  700. x_readonly = x.copy()
  701. x_readonly.flags.writeable = False
  702. peaks, _ = find_peaks(x)
  703. peaks_readonly, _ = find_peaks(x_readonly, **kwargs)
  704. xp_assert_close(peaks, peaks_readonly)
  705. class TestFindPeaksCwt:
  706. def test_find_peaks_exact(self):
  707. """
  708. Generate a series of gaussians and attempt to find the peak locations.
  709. """
  710. sigmas = [5.0, 3.0, 10.0, 20.0, 10.0, 50.0]
  711. num_points = 500
  712. test_data, act_locs = _gen_gaussians_even(sigmas, num_points)
  713. widths = np.arange(0.1, max(sigmas))
  714. found_locs = find_peaks_cwt(test_data, widths, gap_thresh=2, min_snr=0,
  715. min_length=None)
  716. xp_assert_equal(found_locs, act_locs,
  717. check_dtype=False,
  718. err_msg="Found maximum locations did not equal those expected"
  719. )
  720. def test_find_peaks_withnoise(self):
  721. """
  722. Verify that peak locations are (approximately) found
  723. for a series of gaussians with added noise.
  724. """
  725. sigmas = [5.0, 3.0, 10.0, 20.0, 10.0, 50.0]
  726. num_points = 500
  727. test_data, act_locs = _gen_gaussians_even(sigmas, num_points)
  728. widths = np.arange(0.1, max(sigmas))
  729. noise_amp = 0.07
  730. rng = np.random.default_rng(18181911)
  731. test_data += (rng.random(num_points) - 0.5)*(2*noise_amp)
  732. found_locs = find_peaks_cwt(test_data, widths, min_length=15,
  733. gap_thresh=1, min_snr=noise_amp / 5)
  734. err_msg ='Different number of peaks found than expected'
  735. assert len(found_locs) == len(act_locs), err_msg
  736. diffs = np.abs(found_locs - act_locs)
  737. max_diffs = np.array(sigmas) / 5
  738. np.testing.assert_array_less(diffs, max_diffs, 'Maximum location differed' +
  739. f'by more than {max_diffs}')
  740. def test_find_peaks_nopeak(self):
  741. """
  742. Verify that no peak is found in
  743. data that's just noise.
  744. """
  745. noise_amp = 1.0
  746. num_points = 100
  747. rng = np.random.RandomState(181819141)
  748. test_data = (rng.rand(num_points) - 0.5)*(2*noise_amp)
  749. widths = np.arange(10, 50)
  750. found_locs = find_peaks_cwt(test_data, widths, min_snr=5, noise_perc=30)
  751. assert len(found_locs) == 0
  752. def test_find_peaks_with_non_default_wavelets(self):
  753. x = gaussian(200, 2)
  754. widths = np.array([1, 2, 3, 4])
  755. a = find_peaks_cwt(x, widths, wavelet=gaussian)
  756. xp_assert_equal(a, np.asarray([100]), check_dtype=False)
  757. def test_find_peaks_window_size(self):
  758. """
  759. Verify that window_size is passed correctly to private function and
  760. affects the result.
  761. """
  762. sigmas = [2.0, 2.0]
  763. num_points = 1000
  764. test_data, act_locs = _gen_gaussians_even(sigmas, num_points)
  765. widths = np.arange(0.1, max(sigmas), 0.2)
  766. noise_amp = 0.05
  767. rng = np.random.RandomState(18181911)
  768. test_data += (rng.rand(num_points) - 0.5)*(2*noise_amp)
  769. # Possibly contrived negative region to throw off peak finding
  770. # when window_size is too large
  771. test_data[250:320] -= 1
  772. found_locs = find_peaks_cwt(test_data, widths, gap_thresh=2, min_snr=3,
  773. min_length=None, window_size=None)
  774. with pytest.raises(AssertionError):
  775. assert found_locs.size == act_locs.size
  776. found_locs = find_peaks_cwt(test_data, widths, gap_thresh=2, min_snr=3,
  777. min_length=None, window_size=20)
  778. assert found_locs.size == act_locs.size
  779. def test_find_peaks_with_one_width(self):
  780. """
  781. Verify that the `width` argument
  782. in `find_peaks_cwt` can be a float
  783. """
  784. xs = np.arange(0, np.pi, 0.05)
  785. test_data = np.sin(xs)
  786. widths = 1
  787. found_locs = find_peaks_cwt(test_data, widths)
  788. np.testing.assert_equal(found_locs, 32)