test_peak_finding.py 33 KB

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