test_peak_finding.py 31 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824
  1. from __future__ import division, print_function, absolute_import
  2. import copy
  3. import numpy as np
  4. from numpy.testing import (
  5. assert_,
  6. assert_equal,
  7. assert_allclose,
  8. assert_array_equal
  9. )
  10. import pytest
  11. from pytest import raises, warns
  12. from scipy._lib.six import xrange
  13. from scipy.signal._peak_finding import (
  14. argrelmax,
  15. argrelmin,
  16. peak_prominences,
  17. peak_widths,
  18. _unpack_condition_args,
  19. find_peaks,
  20. find_peaks_cwt,
  21. _identify_ridge_lines
  22. )
  23. from scipy.signal._peak_finding_utils import _local_maxima_1d, PeakPropertyWarning
  24. def _gen_gaussians(center_locs, sigmas, total_length):
  25. xdata = np.arange(0, total_length).astype(float)
  26. out_data = np.zeros(total_length, dtype=float)
  27. for ind, sigma in enumerate(sigmas):
  28. tmp = (xdata - center_locs[ind]) / sigma
  29. out_data += np.exp(-(tmp**2))
  30. return out_data
  31. def _gen_gaussians_even(sigmas, total_length):
  32. num_peaks = len(sigmas)
  33. delta = total_length / (num_peaks + 1)
  34. center_locs = np.linspace(delta, total_length - delta, num=num_peaks).astype(int)
  35. out_data = _gen_gaussians(center_locs, sigmas, total_length)
  36. return out_data, center_locs
  37. def _gen_ridge_line(start_locs, max_locs, length, distances, gaps):
  38. """
  39. Generate coordinates for a ridge line.
  40. Will be a series of coordinates, starting a start_loc (length 2).
  41. The maximum distance between any adjacent columns will be
  42. `max_distance`, the max distance between adjacent rows
  43. will be `map_gap'.
  44. `max_locs` should be the size of the intended matrix. The
  45. ending coordinates are guaranteed to be less than `max_locs`,
  46. although they may not approach `max_locs` at all.
  47. """
  48. def keep_bounds(num, max_val):
  49. out = max(num, 0)
  50. out = min(out, max_val)
  51. return out
  52. gaps = copy.deepcopy(gaps)
  53. distances = copy.deepcopy(distances)
  54. locs = np.zeros([length, 2], dtype=int)
  55. locs[0, :] = start_locs
  56. total_length = max_locs[0] - start_locs[0] - sum(gaps)
  57. if total_length < length:
  58. raise ValueError('Cannot generate ridge line according to constraints')
  59. dist_int = length / len(distances) - 1
  60. gap_int = length / len(gaps) - 1
  61. for ind in xrange(1, length):
  62. nextcol = locs[ind - 1, 1]
  63. nextrow = locs[ind - 1, 0] + 1
  64. if (ind % dist_int == 0) and (len(distances) > 0):
  65. nextcol += ((-1)**ind)*distances.pop()
  66. if (ind % gap_int == 0) and (len(gaps) > 0):
  67. nextrow += gaps.pop()
  68. nextrow = keep_bounds(nextrow, max_locs[0])
  69. nextcol = keep_bounds(nextcol, max_locs[1])
  70. locs[ind, :] = [nextrow, nextcol]
  71. return [locs[:, 0], locs[:, 1]]
  72. class TestLocalMaxima1d(object):
  73. def test_empty(self):
  74. """Test with empty signal."""
  75. x = np.array([], dtype=np.float64)
  76. for array in _local_maxima_1d(x):
  77. assert_equal(array, np.array([]))
  78. assert_(array.base is None)
  79. def test_linear(self):
  80. """Test with linear signal."""
  81. x = np.linspace(0, 100)
  82. for array in _local_maxima_1d(x):
  83. assert_equal(array, np.array([]))
  84. assert_(array.base is None)
  85. def test_simple(self):
  86. """Test with simple signal."""
  87. x = np.linspace(-10, 10, 50)
  88. x[2::3] += 1
  89. expected = np.arange(2, 50, 3)
  90. for array in _local_maxima_1d(x):
  91. # For plateaus of size 1, the edges are identical with the
  92. # midpoints
  93. assert_equal(array, expected)
  94. assert_(array.base is None)
  95. def test_flat_maxima(self):
  96. """Test if flat maxima are detected correctly."""
  97. x = np.array([-1.3, 0, 1, 0, 2, 2, 0, 3, 3, 3, 2.99, 4, 4, 4, 4, -10,
  98. -5, -5, -5, -5, -5, -10])
  99. midpoints, left_edges, right_edges = _local_maxima_1d(x)
  100. assert_equal(midpoints, np.array([2, 4, 8, 12, 18]))
  101. assert_equal(left_edges, np.array([2, 4, 7, 11, 16]))
  102. assert_equal(right_edges, np.array([2, 5, 9, 14, 20]))
  103. @pytest.mark.parametrize('x', [
  104. np.array([1., 0, 2]),
  105. np.array([3., 3, 0, 4, 4]),
  106. np.array([5., 5, 5, 0, 6, 6, 6]),
  107. ])
  108. def test_signal_edges(self, x):
  109. """Test if behavior on signal edges is correct."""
  110. for array in _local_maxima_1d(x):
  111. assert_equal(array, np.array([]))
  112. assert_(array.base is None)
  113. def test_exceptions(self):
  114. """Test input validation and raised exceptions."""
  115. with raises(ValueError, match="wrong number of dimensions"):
  116. _local_maxima_1d(np.ones((1, 1)))
  117. with raises(ValueError, match="expected 'float64_t'"):
  118. _local_maxima_1d(np.ones(1, dtype=int))
  119. with raises(TypeError, match="list"):
  120. _local_maxima_1d([1., 2.])
  121. with raises(TypeError, match="'x' must not be None"):
  122. _local_maxima_1d(None)
  123. class TestRidgeLines(object):
  124. def test_empty(self):
  125. test_matr = np.zeros([20, 100])
  126. lines = _identify_ridge_lines(test_matr, 2*np.ones(20), 1)
  127. assert_(len(lines) == 0)
  128. def test_minimal(self):
  129. test_matr = np.zeros([20, 100])
  130. test_matr[0, 10] = 1
  131. lines = _identify_ridge_lines(test_matr, 2*np.ones(20), 1)
  132. assert_(len(lines) == 1)
  133. test_matr = np.zeros([20, 100])
  134. test_matr[0:2, 10] = 1
  135. lines = _identify_ridge_lines(test_matr, 2*np.ones(20), 1)
  136. assert_(len(lines) == 1)
  137. def test_single_pass(self):
  138. distances = [0, 1, 2, 5]
  139. gaps = [0, 1, 2, 0, 1]
  140. test_matr = np.zeros([20, 50]) + 1e-12
  141. length = 12
  142. line = _gen_ridge_line([0, 25], test_matr.shape, length, distances, gaps)
  143. test_matr[line[0], line[1]] = 1
  144. max_distances = max(distances)*np.ones(20)
  145. identified_lines = _identify_ridge_lines(test_matr, max_distances, max(gaps) + 1)
  146. assert_array_equal(identified_lines, [line])
  147. def test_single_bigdist(self):
  148. distances = [0, 1, 2, 5]
  149. gaps = [0, 1, 2, 4]
  150. test_matr = np.zeros([20, 50])
  151. length = 12
  152. line = _gen_ridge_line([0, 25], test_matr.shape, length, distances, gaps)
  153. test_matr[line[0], line[1]] = 1
  154. max_dist = 3
  155. max_distances = max_dist*np.ones(20)
  156. #This should get 2 lines, since the distance is too large
  157. identified_lines = _identify_ridge_lines(test_matr, max_distances, 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 = max_dist*np.ones(20)
  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 = max_dist*np.ones(50)
  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(object):
  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. assert_equal(len(i), 1)
  210. assert_array_equal(i[0], empty_array)
  211. z2 = np.zeros((3,5))
  212. row, col = argrelmin(z2, axis=0)
  213. assert_array_equal(row, empty_array)
  214. assert_array_equal(col, empty_array)
  215. row, col = argrelmin(z2, axis=1)
  216. assert_array_equal(row, empty_array)
  217. assert_array_equal(col, empty_array)
  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. assert_equal(row[order], [1, 2, 3])
  230. assert_equal(col[order], [4, 0, 1])
  231. row, col = argrelmax(x, axis=1)
  232. order = np.argsort(row)
  233. assert_equal(row[order], [0, 3, 4])
  234. assert_equal(col[order], [3, 1, 2])
  235. row, col = argrelmin(x, axis=0)
  236. order = np.argsort(row)
  237. assert_equal(row[order], [1, 2, 3])
  238. assert_equal(col[order], [1, 2, 3])
  239. row, col = argrelmin(x, axis=1)
  240. order = np.argsort(row)
  241. assert_equal(row[order], [1, 2, 3])
  242. assert_equal(col[order], [1, 2, 3])
  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 xrange(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(object):
  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. assert_equal(out[0], proms)
  290. assert_equal(out[1], lbases)
  291. assert_equal(out[2], rbases)
  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. assert_equal(proms, [2, 2, 2])
  301. assert_equal(lbases, [0, 0, 0])
  302. assert_equal(rbases, [6, 6, 6])
  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. assert_equal(proms, [1, 1, 1])
  308. assert_equal(lbases, peaks - 1)
  309. assert_equal(rbases, peaks + 1)
  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. assert_equal(proms, [9, 9, 2])
  318. assert_equal(lbases, [0, 0, 3])
  319. assert_equal(rbases, [3, 3, 5])
  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. assert_equal(peak_prominences(x, peak), [3., 0, 6])
  328. for wlen, i in [(8, 0), (7, 0), (6, 0), (5, 1), (3.2, 1), (3, 2), (1.1, 2)]:
  329. assert_equal(peak_prominences(x, peak, wlen), [3. - i, 0 + i, 6 - i])
  330. def test_exceptions(self):
  331. """
  332. Verify that exceptions and warnings are raised.
  333. """
  334. # x with dimension > 1
  335. with raises(ValueError, match='1D array'):
  336. peak_prominences([[0, 1, 1, 0]], [1, 2])
  337. # peaks with dimension > 1
  338. with raises(ValueError, match='1D array'):
  339. peak_prominences([0, 1, 1, 0], [[1, 2]])
  340. # x with dimension < 1
  341. with raises(ValueError, match='1D array'):
  342. peak_prominences(3, [0,])
  343. # empty x with supplied
  344. with raises(ValueError, match='not a valid index'):
  345. peak_prominences([], [0])
  346. # invalid indices with non-empty x
  347. for p in [-100, -1, 3, 1000]:
  348. with raises(ValueError, match='not a valid index'):
  349. peak_prominences([1, 0, 2], [p])
  350. # peaks is not cast-able to np.intp
  351. with raises(TypeError, match='cannot safely cast'):
  352. peak_prominences([0, 1, 1, 0], [1.1, 2.3])
  353. # wlen < 3
  354. with raises(ValueError, match='wlen'):
  355. peak_prominences(np.arange(10), [3, 5], wlen=1)
  356. def test_warnings(self):
  357. """
  358. Verify that appropriate warnings are raised.
  359. """
  360. msg = "some peaks have a prominence of 0"
  361. for p in [0, 1, 2]:
  362. with warns(PeakPropertyWarning, match=msg):
  363. peak_prominences([1, 0, 2], [p,])
  364. with warns(PeakPropertyWarning, match=msg):
  365. peak_prominences([0, 1, 1, 1, 0], [2], wlen=2)
  366. class TestPeakWidths(object):
  367. def test_empty(self):
  368. """
  369. Test if an empty array is returned if no peaks are provided.
  370. """
  371. widths = peak_widths([], [])[0]
  372. assert_(isinstance(widths, np.ndarray))
  373. assert_equal(widths.size, 0)
  374. widths = peak_widths([1, 2, 3], [])[0]
  375. assert_(isinstance(widths, np.ndarray))
  376. assert_equal(widths.size, 0)
  377. out = peak_widths([], [])
  378. for arr in out:
  379. assert_(isinstance(arr, np.ndarray))
  380. assert_equal(arr.size, 0)
  381. @pytest.mark.filterwarnings("ignore:some peaks have a width of 0")
  382. def test_basic(self):
  383. """
  384. Test a simple use case with easy to verify results at different relative
  385. heights.
  386. """
  387. x = np.array([1, 0, 1, 2, 1, 0, -1])
  388. prominence = 2
  389. for rel_height, width_true, lip_true, rip_true in [
  390. (0., 0., 3., 3.), # raises warning
  391. (0.25, 1., 2.5, 3.5),
  392. (0.5, 2., 2., 4.),
  393. (0.75, 3., 1.5, 4.5),
  394. (1., 4., 1., 5.),
  395. (2., 5., 1., 6.),
  396. (3., 5., 1., 6.)
  397. ]:
  398. width_calc, height, lip_calc, rip_calc = peak_widths(
  399. x, [3], rel_height)
  400. assert_allclose(width_calc, width_true)
  401. assert_allclose(height, 2 - rel_height * prominence)
  402. assert_allclose(lip_calc, lip_true)
  403. assert_allclose(rip_calc, rip_true)
  404. def test_non_contiguous(self):
  405. """
  406. Test with non-C-contiguous input arrays.
  407. """
  408. x = np.repeat([0, 100, 50], 4)
  409. peaks = np.repeat([1], 3)
  410. result = peak_widths(x[::4], peaks[::3])
  411. assert_equal(result, [0.75, 75, 0.75, 1.5])
  412. def test_exceptions(self):
  413. """
  414. Verify that argument validation works as intended.
  415. """
  416. with raises(ValueError, match='1D array'):
  417. # x with dimension > 1
  418. peak_widths(np.zeros((3, 4)), np.ones(3))
  419. with raises(ValueError, match='1D array'):
  420. # x with dimension < 1
  421. peak_widths(3, [0])
  422. with raises(ValueError, match='1D array'):
  423. # peaks with dimension > 1
  424. peak_widths(np.arange(10), np.ones((3, 2), dtype=np.intp))
  425. with raises(ValueError, match='1D array'):
  426. # peaks with dimension < 1
  427. peak_widths(np.arange(10), 3)
  428. with raises(ValueError, match='not a valid index'):
  429. # peak pos exceeds x.size
  430. peak_widths(np.arange(10), [8, 11])
  431. with raises(ValueError, match='not a valid index'):
  432. # empty x with peaks supplied
  433. peak_widths([], [1, 2])
  434. with raises(TypeError, match='cannot safely cast'):
  435. # peak cannot be safely casted to intp
  436. peak_widths(np.arange(10), [1.1, 2.3])
  437. with raises(ValueError, match='rel_height'):
  438. # rel_height is < 0
  439. peak_widths([0, 1, 0, 1, 0], [1, 3], rel_height=-1)
  440. with raises(TypeError, match='None'):
  441. # prominence data contains None
  442. peak_widths([1, 2, 1], [1], prominence_data=(None, None, None))
  443. def test_warnings(self):
  444. """
  445. Verify that appropriate warnings are raised.
  446. """
  447. msg = "some peaks have a width of 0"
  448. with warns(PeakPropertyWarning, match=msg):
  449. # Case: rel_height is 0
  450. peak_widths([0, 1, 0], [1], rel_height=0)
  451. with warns(PeakPropertyWarning, match=msg):
  452. # Case: prominence is 0 and bases are identical
  453. peak_widths(
  454. [0, 1, 1, 1, 0], [2],
  455. prominence_data=(np.array([0.], np.float64),
  456. np.array([2], np.intp),
  457. np.array([2], np.intp))
  458. )
  459. def test_mismatching_prominence_data(self):
  460. """Test with mismatching peak and / or prominence data."""
  461. x = [0, 1, 0]
  462. peak = [1]
  463. for i, (prominences, left_bases, right_bases) in enumerate([
  464. ((1.,), (-1,), (2,)), # left base not in x
  465. ((1.,), (0,), (3,)), # right base not in x
  466. ((1.,), (2,), (0,)), # swapped bases same as peak
  467. ((1., 1.), (0, 0), (2, 2)), # array shapes don't match peaks
  468. ((1., 1.), (0,), (2,)), # arrays with different shapes
  469. ((1.,), (0, 0), (2,)), # arrays with different shapes
  470. ((1.,), (0,), (2, 2)) # arrays with different shapes
  471. ]):
  472. # Make sure input is matches output of signal.peak_prominences
  473. prominence_data = (np.array(prominences, dtype=np.float64),
  474. np.array(left_bases, dtype=np.intp),
  475. np.array(right_bases, dtype=np.intp))
  476. # Test for correct exception
  477. if i < 3:
  478. match = "prominence data is invalid for peak"
  479. else:
  480. match = "arrays in `prominence_data` must have the same shape"
  481. with raises(ValueError, match=match):
  482. peak_widths(x, peak, prominence_data=prominence_data)
  483. @pytest.mark.filterwarnings("ignore:some peaks have a width of 0")
  484. def test_intersection_rules(self):
  485. """Test if x == eval_height counts as an intersection."""
  486. # Flatt peak with two possible intersection points if evaluated at 1
  487. x = [0, 1, 2, 1, 3, 3, 3, 1, 2, 1, 0]
  488. # relative height is 0 -> width is 0 as well, raises warning
  489. assert_allclose(peak_widths(x, peaks=[5], rel_height=0),
  490. [(0.,), (3.,), (5.,), (5.,)])
  491. # width_height == x counts as intersection -> nearest 1 is chosen
  492. assert_allclose(peak_widths(x, peaks=[5], rel_height=2/3),
  493. [(4.,), (1.,), (3.,), (7.,)])
  494. def test_unpack_condition_args():
  495. """
  496. Verify parsing of condition arguments for `scipy.signal.find_peaks` function.
  497. """
  498. x = np.arange(10)
  499. amin_true = x
  500. amax_true = amin_true + 10
  501. peaks = amin_true[1::2]
  502. # Test unpacking with None or interval
  503. assert_((None, None) == _unpack_condition_args((None, None), x, peaks))
  504. assert_((1, None) == _unpack_condition_args(1, x, peaks))
  505. assert_((1, None) == _unpack_condition_args((1, None), x, peaks))
  506. assert_((None, 2) == _unpack_condition_args((None, 2), x, peaks))
  507. assert_((3., 4.5) == _unpack_condition_args((3., 4.5), x, peaks))
  508. # Test if borders are correctly reduced with `peaks`
  509. amin_calc, amax_calc = _unpack_condition_args((amin_true, amax_true), x, peaks)
  510. assert_equal(amin_calc, amin_true[peaks])
  511. assert_equal(amax_calc, amax_true[peaks])
  512. # Test raises if array borders don't match x
  513. with raises(ValueError, match="array size of lower"):
  514. _unpack_condition_args(amin_true, np.arange(11), peaks)
  515. with raises(ValueError, match="array size of upper"):
  516. _unpack_condition_args((None, amin_true), np.arange(11), peaks)
  517. class TestFindPeaks(object):
  518. # Keys of optionally returned properties
  519. property_keys = {'peak_heights', 'left_thresholds', 'right_thresholds',
  520. 'prominences', 'left_bases', 'right_bases', 'widths',
  521. 'width_heights', 'left_ips', 'right_ips'}
  522. def test_constant(self):
  523. """
  524. Test behavior for signal without local maxima.
  525. """
  526. open_interval = (None, None)
  527. peaks, props = find_peaks(np.ones(10),
  528. height=open_interval, threshold=open_interval,
  529. prominence=open_interval, width=open_interval)
  530. assert_(peaks.size == 0)
  531. for key in self.property_keys:
  532. assert_(props[key].size == 0)
  533. def test_plateau_size(self):
  534. """
  535. Test plateau size condition for peaks.
  536. """
  537. # Prepare signal with peaks with peak_height == plateau_size
  538. plateau_sizes = np.array([1, 2, 3, 4, 8, 20, 111])
  539. x = np.zeros(plateau_sizes.size * 2 + 1)
  540. x[1::2] = plateau_sizes
  541. repeats = np.ones(x.size, dtype=int)
  542. repeats[1::2] = x[1::2]
  543. x = np.repeat(x, repeats)
  544. # Test full output
  545. peaks, props = find_peaks(x, plateau_size=(None, None))
  546. assert_equal(peaks, [1, 3, 7, 11, 18, 33, 100])
  547. assert_equal(props["plateau_sizes"], plateau_sizes)
  548. assert_equal(props["left_edges"], peaks - (plateau_sizes - 1) // 2)
  549. assert_equal(props["right_edges"], peaks + plateau_sizes // 2)
  550. # Test conditions
  551. assert_equal(find_peaks(x, plateau_size=4)[0], [11, 18, 33, 100])
  552. assert_equal(find_peaks(x, plateau_size=(None, 3.5))[0], [1, 3, 7])
  553. assert_equal(find_peaks(x, plateau_size=(5, 50))[0], [18, 33])
  554. def test_height_condition(self):
  555. """
  556. Test height condition for peaks.
  557. """
  558. x = (0., 1/3, 0., 2.5, 0, 4., 0)
  559. peaks, props = find_peaks(x, height=(None, None))
  560. assert_equal(peaks, np.array([1, 3, 5]))
  561. assert_equal(props['peak_heights'], np.array([1/3, 2.5, 4.]))
  562. assert_equal(find_peaks(x, height=0.5)[0], np.array([3, 5]))
  563. assert_equal(find_peaks(x, height=(None, 3))[0], np.array([1, 3]))
  564. assert_equal(find_peaks(x, height=(2, 3))[0], np.array([3]))
  565. def test_threshold_condition(self):
  566. """
  567. Test threshold condition for peaks.
  568. """
  569. x = (0, 2, 1, 4, -1)
  570. peaks, props = find_peaks(x, threshold=(None, None))
  571. assert_equal(peaks, np.array([1, 3]))
  572. assert_equal(props['left_thresholds'], np.array([2, 3]))
  573. assert_equal(props['right_thresholds'], np.array([1, 5]))
  574. assert_equal(find_peaks(x, threshold=2)[0], np.array([3]))
  575. assert_equal(find_peaks(x, threshold=3.5)[0], np.array([]))
  576. assert_equal(find_peaks(x, threshold=(None, 5))[0], np.array([1, 3]))
  577. assert_equal(find_peaks(x, threshold=(None, 4))[0], np.array([1]))
  578. assert_equal(find_peaks(x, threshold=(2, 4))[0], np.array([]))
  579. def test_distance_condition(self):
  580. """
  581. Test distance condition for peaks.
  582. """
  583. # Peaks of different height with constant distance 3
  584. peaks_all = np.arange(1, 21, 3)
  585. x = np.zeros(21)
  586. x[peaks_all] += np.linspace(1, 2, peaks_all.size)
  587. # Test if peaks with "minimal" distance are still selected (distance = 3)
  588. assert_equal(find_peaks(x, distance=3)[0], peaks_all)
  589. # Select every second peak (distance > 3)
  590. peaks_subset = find_peaks(x, distance=3.0001)[0]
  591. # Test if peaks_subset is subset of peaks_all
  592. assert_(
  593. np.setdiff1d(peaks_subset, peaks_all, assume_unique=True).size == 0
  594. )
  595. # Test if every second peak was removed
  596. assert_equal(np.diff(peaks_subset), 6)
  597. # Test priority of peak removal
  598. x = [-2, 1, -1, 0, -3]
  599. peaks_subset = find_peaks(x, distance=10)[0] # use distance > x size
  600. assert_(peaks_subset.size == 1 and peaks_subset[0] == 1)
  601. def test_prominence_condition(self):
  602. """
  603. Test prominence condition for peaks.
  604. """
  605. x = np.linspace(0, 10, 100)
  606. peaks_true = np.arange(1, 99, 2)
  607. offset = np.linspace(1, 10, peaks_true.size)
  608. x[peaks_true] += offset
  609. prominences = x[peaks_true] - x[peaks_true + 1]
  610. interval = (3, 9)
  611. keep = np.nonzero(
  612. (interval[0] <= prominences) & (prominences <= interval[1]))
  613. peaks_calc, properties = find_peaks(x, prominence=interval)
  614. assert_equal(peaks_calc, peaks_true[keep])
  615. assert_equal(properties['prominences'], prominences[keep])
  616. assert_equal(properties['left_bases'], 0)
  617. assert_equal(properties['right_bases'], peaks_true[keep] + 1)
  618. def test_width_condition(self):
  619. """
  620. Test width condition for peaks.
  621. """
  622. x = np.array([1, 0, 1, 2, 1, 0, -1, 4, 0])
  623. peaks, props = find_peaks(x, width=(None, 2), rel_height=0.75)
  624. assert_equal(peaks.size, 1)
  625. assert_equal(peaks, 7)
  626. assert_allclose(props['widths'], 1.35)
  627. assert_allclose(props['width_heights'], 1.)
  628. assert_allclose(props['left_ips'], 6.4)
  629. assert_allclose(props['right_ips'], 7.75)
  630. def test_properties(self):
  631. """
  632. Test returned properties.
  633. """
  634. open_interval = (None, None)
  635. x = [0, 1, 0, 2, 1.5, 0, 3, 0, 5, 9]
  636. peaks, props = find_peaks(x,
  637. height=open_interval, threshold=open_interval,
  638. prominence=open_interval, width=open_interval)
  639. assert_(len(props) == len(self.property_keys))
  640. for key in self.property_keys:
  641. assert_(peaks.size == props[key].size)
  642. def test_raises(self):
  643. """
  644. Test exceptions raised by function.
  645. """
  646. with raises(ValueError, match="1D array"):
  647. find_peaks(np.array(1))
  648. with raises(ValueError, match="1D array"):
  649. find_peaks(np.ones((2, 2)))
  650. with raises(ValueError, match="distance"):
  651. find_peaks(np.arange(10), distance=-1)
  652. @pytest.mark.filterwarnings("ignore:some peaks have a prominence of 0",
  653. "ignore:some peaks have a width of 0")
  654. def test_wlen_smaller_plateau(self):
  655. """
  656. Test behavior of prominence and width calculation if the given window
  657. length is smaller than a peak's plateau size.
  658. Regression test for gh-9110.
  659. """
  660. peaks, props = find_peaks([0, 1, 1, 1, 0], prominence=(None, None),
  661. width=(None, None), wlen=2)
  662. assert_equal(peaks, 2)
  663. assert_equal(props["prominences"], 0)
  664. assert_equal(props["widths"], 0)
  665. assert_equal(props["width_heights"], 1)
  666. for key in ("left_bases", "right_bases", "left_ips", "right_ips"):
  667. assert_equal(props[key], peaks)
  668. class TestFindPeaksCwt(object):
  669. def test_find_peaks_exact(self):
  670. """
  671. Generate a series of gaussians and attempt to find the peak locations.
  672. """
  673. sigmas = [5.0, 3.0, 10.0, 20.0, 10.0, 50.0]
  674. num_points = 500
  675. test_data, act_locs = _gen_gaussians_even(sigmas, num_points)
  676. widths = np.arange(0.1, max(sigmas))
  677. found_locs = find_peaks_cwt(test_data, widths, gap_thresh=2, min_snr=0,
  678. min_length=None)
  679. np.testing.assert_array_equal(found_locs, act_locs,
  680. "Found maximum locations did not equal those expected")
  681. def test_find_peaks_withnoise(self):
  682. """
  683. Verify that peak locations are (approximately) found
  684. for a series of gaussians with added noise.
  685. """
  686. sigmas = [5.0, 3.0, 10.0, 20.0, 10.0, 50.0]
  687. num_points = 500
  688. test_data, act_locs = _gen_gaussians_even(sigmas, num_points)
  689. widths = np.arange(0.1, max(sigmas))
  690. noise_amp = 0.07
  691. np.random.seed(18181911)
  692. test_data += (np.random.rand(num_points) - 0.5)*(2*noise_amp)
  693. found_locs = find_peaks_cwt(test_data, widths, min_length=15,
  694. gap_thresh=1, min_snr=noise_amp / 5)
  695. np.testing.assert_equal(len(found_locs), len(act_locs), 'Different number' +
  696. 'of peaks found than expected')
  697. diffs = np.abs(found_locs - act_locs)
  698. max_diffs = np.array(sigmas) / 5
  699. np.testing.assert_array_less(diffs, max_diffs, 'Maximum location differed' +
  700. 'by more than %s' % (max_diffs))
  701. def test_find_peaks_nopeak(self):
  702. """
  703. Verify that no peak is found in
  704. data that's just noise.
  705. """
  706. noise_amp = 1.0
  707. num_points = 100
  708. np.random.seed(181819141)
  709. test_data = (np.random.rand(num_points) - 0.5)*(2*noise_amp)
  710. widths = np.arange(10, 50)
  711. found_locs = find_peaks_cwt(test_data, widths, min_snr=5, noise_perc=30)
  712. np.testing.assert_equal(len(found_locs), 0)