test_vq.py 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336
  1. import warnings
  2. import sys
  3. import numpy as np
  4. from numpy.testing import (assert_array_equal, assert_array_almost_equal,
  5. assert_allclose, assert_equal, assert_,
  6. suppress_warnings)
  7. import pytest
  8. from pytest import raises as assert_raises
  9. from scipy.cluster.vq import (kmeans, kmeans2, py_vq, vq, whiten,
  10. ClusterError, _krandinit)
  11. from scipy.cluster import _vq
  12. from scipy.sparse._sputils import matrix
  13. TESTDATA_2D = np.array([
  14. -2.2, 1.17, -1.63, 1.69, -2.04, 4.38, -3.09, 0.95, -1.7, 4.79, -1.68, 0.68,
  15. -2.26, 3.34, -2.29, 2.55, -1.72, -0.72, -1.99, 2.34, -2.75, 3.43, -2.45,
  16. 2.41, -4.26, 3.65, -1.57, 1.87, -1.96, 4.03, -3.01, 3.86, -2.53, 1.28,
  17. -4.0, 3.95, -1.62, 1.25, -3.42, 3.17, -1.17, 0.12, -3.03, -0.27, -2.07,
  18. -0.55, -1.17, 1.34, -2.82, 3.08, -2.44, 0.24, -1.71, 2.48, -5.23, 4.29,
  19. -2.08, 3.69, -1.89, 3.62, -2.09, 0.26, -0.92, 1.07, -2.25, 0.88, -2.25,
  20. 2.02, -4.31, 3.86, -2.03, 3.42, -2.76, 0.3, -2.48, -0.29, -3.42, 3.21,
  21. -2.3, 1.73, -2.84, 0.69, -1.81, 2.48, -5.24, 4.52, -2.8, 1.31, -1.67,
  22. -2.34, -1.18, 2.17, -2.17, 2.82, -1.85, 2.25, -2.45, 1.86, -6.79, 3.94,
  23. -2.33, 1.89, -1.55, 2.08, -1.36, 0.93, -2.51, 2.74, -2.39, 3.92, -3.33,
  24. 2.99, -2.06, -0.9, -2.83, 3.35, -2.59, 3.05, -2.36, 1.85, -1.69, 1.8,
  25. -1.39, 0.66, -2.06, 0.38, -1.47, 0.44, -4.68, 3.77, -5.58, 3.44, -2.29,
  26. 2.24, -1.04, -0.38, -1.85, 4.23, -2.88, 0.73, -2.59, 1.39, -1.34, 1.75,
  27. -1.95, 1.3, -2.45, 3.09, -1.99, 3.41, -5.55, 5.21, -1.73, 2.52, -2.17,
  28. 0.85, -2.06, 0.49, -2.54, 2.07, -2.03, 1.3, -3.23, 3.09, -1.55, 1.44,
  29. -0.81, 1.1, -2.99, 2.92, -1.59, 2.18, -2.45, -0.73, -3.12, -1.3, -2.83,
  30. 0.2, -2.77, 3.24, -1.98, 1.6, -4.59, 3.39, -4.85, 3.75, -2.25, 1.71, -3.28,
  31. 3.38, -1.74, 0.88, -2.41, 1.92, -2.24, 1.19, -2.48, 1.06, -1.68, -0.62,
  32. -1.3, 0.39, -1.78, 2.35, -3.54, 2.44, -1.32, 0.66, -2.38, 2.76, -2.35,
  33. 3.95, -1.86, 4.32, -2.01, -1.23, -1.79, 2.76, -2.13, -0.13, -5.25, 3.84,
  34. -2.24, 1.59, -4.85, 2.96, -2.41, 0.01, -0.43, 0.13, -3.92, 2.91, -1.75,
  35. -0.53, -1.69, 1.69, -1.09, 0.15, -2.11, 2.17, -1.53, 1.22, -2.1, -0.86,
  36. -2.56, 2.28, -3.02, 3.33, -1.12, 3.86, -2.18, -1.19, -3.03, 0.79, -0.83,
  37. 0.97, -3.19, 1.45, -1.34, 1.28, -2.52, 4.22, -4.53, 3.22, -1.97, 1.75,
  38. -2.36, 3.19, -0.83, 1.53, -1.59, 1.86, -2.17, 2.3, -1.63, 2.71, -2.03,
  39. 3.75, -2.57, -0.6, -1.47, 1.33, -1.95, 0.7, -1.65, 1.27, -1.42, 1.09, -3.0,
  40. 3.87, -2.51, 3.06, -2.6, 0.74, -1.08, -0.03, -2.44, 1.31, -2.65, 2.99,
  41. -1.84, 1.65, -4.76, 3.75, -2.07, 3.98, -2.4, 2.67, -2.21, 1.49, -1.21,
  42. 1.22, -5.29, 2.38, -2.85, 2.28, -5.6, 3.78, -2.7, 0.8, -1.81, 3.5, -3.75,
  43. 4.17, -1.29, 2.99, -5.92, 3.43, -1.83, 1.23, -1.24, -1.04, -2.56, 2.37,
  44. -3.26, 0.39, -4.63, 2.51, -4.52, 3.04, -1.7, 0.36, -1.41, 0.04, -2.1, 1.0,
  45. -1.87, 3.78, -4.32, 3.59, -2.24, 1.38, -1.99, -0.22, -1.87, 1.95, -0.84,
  46. 2.17, -5.38, 3.56, -1.27, 2.9, -1.79, 3.31, -5.47, 3.85, -1.44, 3.69,
  47. -2.02, 0.37, -1.29, 0.33, -2.34, 2.56, -1.74, -1.27, -1.97, 1.22, -2.51,
  48. -0.16, -1.64, -0.96, -2.99, 1.4, -1.53, 3.31, -2.24, 0.45, -2.46, 1.71,
  49. -2.88, 1.56, -1.63, 1.46, -1.41, 0.68, -1.96, 2.76, -1.61,
  50. 2.11]).reshape((200, 2))
  51. # Global data
  52. X = np.array([[3.0, 3], [4, 3], [4, 2],
  53. [9, 2], [5, 1], [6, 2], [9, 4],
  54. [5, 2], [5, 4], [7, 4], [6, 5]])
  55. CODET1 = np.array([[3.0000, 3.0000],
  56. [6.2000, 4.0000],
  57. [5.8000, 1.8000]])
  58. CODET2 = np.array([[11.0/3, 8.0/3],
  59. [6.7500, 4.2500],
  60. [6.2500, 1.7500]])
  61. LABEL1 = np.array([0, 1, 2, 2, 2, 2, 1, 2, 1, 1, 1])
  62. class TestWhiten:
  63. def test_whiten(self):
  64. desired = np.array([[5.08738849, 2.97091878],
  65. [3.19909255, 0.69660580],
  66. [4.51041982, 0.02640918],
  67. [4.38567074, 0.95120889],
  68. [2.32191480, 1.63195503]])
  69. for tp in np.array, matrix:
  70. obs = tp([[0.98744510, 0.82766775],
  71. [0.62093317, 0.19406729],
  72. [0.87545741, 0.00735733],
  73. [0.85124403, 0.26499712],
  74. [0.45067590, 0.45464607]])
  75. assert_allclose(whiten(obs), desired, rtol=1e-5)
  76. def test_whiten_zero_std(self):
  77. desired = np.array([[0., 1.0, 2.86666544],
  78. [0., 1.0, 1.32460034],
  79. [0., 1.0, 3.74382172]])
  80. for tp in np.array, matrix:
  81. obs = tp([[0., 1., 0.74109533],
  82. [0., 1., 0.34243798],
  83. [0., 1., 0.96785929]])
  84. with warnings.catch_warnings(record=True) as w:
  85. warnings.simplefilter('always')
  86. assert_allclose(whiten(obs), desired, rtol=1e-5)
  87. assert_equal(len(w), 1)
  88. assert_(issubclass(w[-1].category, RuntimeWarning))
  89. def test_whiten_not_finite(self):
  90. for tp in np.array, matrix:
  91. for bad_value in np.nan, np.inf, -np.inf:
  92. obs = tp([[0.98744510, bad_value],
  93. [0.62093317, 0.19406729],
  94. [0.87545741, 0.00735733],
  95. [0.85124403, 0.26499712],
  96. [0.45067590, 0.45464607]])
  97. assert_raises(ValueError, whiten, obs)
  98. class TestVq:
  99. def test_py_vq(self):
  100. initc = np.concatenate(([[X[0]], [X[1]], [X[2]]]))
  101. for tp in np.array, matrix:
  102. label1 = py_vq(tp(X), tp(initc))[0]
  103. assert_array_equal(label1, LABEL1)
  104. def test_vq(self):
  105. initc = np.concatenate(([[X[0]], [X[1]], [X[2]]]))
  106. for tp in np.array, matrix:
  107. label1, dist = _vq.vq(tp(X), tp(initc))
  108. assert_array_equal(label1, LABEL1)
  109. tlabel1, tdist = vq(tp(X), tp(initc))
  110. def test_vq_1d(self):
  111. # Test special rank 1 vq algo, python implementation.
  112. data = X[:, 0]
  113. initc = data[:3]
  114. a, b = _vq.vq(data, initc)
  115. ta, tb = py_vq(data[:, np.newaxis], initc[:, np.newaxis])
  116. assert_array_equal(a, ta)
  117. assert_array_equal(b, tb)
  118. def test__vq_sametype(self):
  119. a = np.array([1.0, 2.0], dtype=np.float64)
  120. b = a.astype(np.float32)
  121. assert_raises(TypeError, _vq.vq, a, b)
  122. def test__vq_invalid_type(self):
  123. a = np.array([1, 2], dtype=int)
  124. assert_raises(TypeError, _vq.vq, a, a)
  125. def test_vq_large_nfeat(self):
  126. X = np.random.rand(20, 20)
  127. code_book = np.random.rand(3, 20)
  128. codes0, dis0 = _vq.vq(X, code_book)
  129. codes1, dis1 = py_vq(X, code_book)
  130. assert_allclose(dis0, dis1, 1e-5)
  131. assert_array_equal(codes0, codes1)
  132. X = X.astype(np.float32)
  133. code_book = code_book.astype(np.float32)
  134. codes0, dis0 = _vq.vq(X, code_book)
  135. codes1, dis1 = py_vq(X, code_book)
  136. assert_allclose(dis0, dis1, 1e-5)
  137. assert_array_equal(codes0, codes1)
  138. def test_vq_large_features(self):
  139. X = np.random.rand(10, 5) * 1000000
  140. code_book = np.random.rand(2, 5) * 1000000
  141. codes0, dis0 = _vq.vq(X, code_book)
  142. codes1, dis1 = py_vq(X, code_book)
  143. assert_allclose(dis0, dis1, 1e-5)
  144. assert_array_equal(codes0, codes1)
  145. class TestKMean:
  146. def test_large_features(self):
  147. # Generate a data set with large values, and run kmeans on it to
  148. # (regression for 1077).
  149. d = 300
  150. n = 100
  151. m1 = np.random.randn(d)
  152. m2 = np.random.randn(d)
  153. x = 10000 * np.random.randn(n, d) - 20000 * m1
  154. y = 10000 * np.random.randn(n, d) + 20000 * m2
  155. data = np.empty((x.shape[0] + y.shape[0], d), np.double)
  156. data[:x.shape[0]] = x
  157. data[x.shape[0]:] = y
  158. kmeans(data, 2)
  159. def test_kmeans_simple(self):
  160. np.random.seed(54321)
  161. initc = np.concatenate(([[X[0]], [X[1]], [X[2]]]))
  162. for tp in np.array, matrix:
  163. code1 = kmeans(tp(X), tp(initc), iter=1)[0]
  164. assert_array_almost_equal(code1, CODET2)
  165. def test_kmeans_lost_cluster(self):
  166. # This will cause kmeans to have a cluster with no points.
  167. data = TESTDATA_2D
  168. initk = np.array([[-1.8127404, -0.67128041],
  169. [2.04621601, 0.07401111],
  170. [-2.31149087, -0.05160469]])
  171. kmeans(data, initk)
  172. with suppress_warnings() as sup:
  173. sup.filter(UserWarning,
  174. "One of the clusters is empty. Re-run kmeans with a "
  175. "different initialization")
  176. kmeans2(data, initk, missing='warn')
  177. assert_raises(ClusterError, kmeans2, data, initk, missing='raise')
  178. def test_kmeans2_simple(self):
  179. np.random.seed(12345678)
  180. initc = np.concatenate(([[X[0]], [X[1]], [X[2]]]))
  181. for tp in np.array, matrix:
  182. code1 = kmeans2(tp(X), tp(initc), iter=1)[0]
  183. code2 = kmeans2(tp(X), tp(initc), iter=2)[0]
  184. assert_array_almost_equal(code1, CODET1)
  185. assert_array_almost_equal(code2, CODET2)
  186. def test_kmeans2_rank1(self):
  187. data = TESTDATA_2D
  188. data1 = data[:, 0]
  189. initc = data1[:3]
  190. code = initc.copy()
  191. kmeans2(data1, code, iter=1)[0]
  192. kmeans2(data1, code, iter=2)[0]
  193. def test_kmeans2_rank1_2(self):
  194. data = TESTDATA_2D
  195. data1 = data[:, 0]
  196. kmeans2(data1, 2, iter=1)
  197. def test_kmeans2_high_dim(self):
  198. # test kmeans2 when the number of dimensions exceeds the number
  199. # of input points
  200. data = TESTDATA_2D
  201. data = data.reshape((20, 20))[:10]
  202. kmeans2(data, 2)
  203. def test_kmeans2_init(self):
  204. np.random.seed(12345)
  205. data = TESTDATA_2D
  206. kmeans2(data, 3, minit='points')
  207. kmeans2(data[:, :1], 3, minit='points') # special case (1-D)
  208. kmeans2(data, 3, minit='++')
  209. kmeans2(data[:, :1], 3, minit='++') # special case (1-D)
  210. # minit='random' can give warnings, filter those
  211. with suppress_warnings() as sup:
  212. sup.filter(message="One of the clusters is empty. Re-run.")
  213. kmeans2(data, 3, minit='random')
  214. kmeans2(data[:, :1], 3, minit='random') # special case (1-D)
  215. @pytest.mark.skipif(sys.platform == 'win32',
  216. reason='Fails with MemoryError in Wine.')
  217. def test_krandinit(self):
  218. data = TESTDATA_2D
  219. datas = [data.reshape((200, 2)), data.reshape((20, 20))[:10]]
  220. k = int(1e6)
  221. for data in datas:
  222. # check that np.random.Generator can be used (numpy >= 1.17)
  223. if hasattr(np.random, 'default_rng'):
  224. rng = np.random.default_rng(1234)
  225. else:
  226. rng = np.random.RandomState(1234)
  227. init = _krandinit(data, k, rng)
  228. orig_cov = np.cov(data, rowvar=0)
  229. init_cov = np.cov(init, rowvar=0)
  230. assert_allclose(orig_cov, init_cov, atol=1e-2)
  231. def test_kmeans2_empty(self):
  232. # Regression test for gh-1032.
  233. assert_raises(ValueError, kmeans2, [], 2)
  234. def test_kmeans_0k(self):
  235. # Regression test for gh-1073: fail when k arg is 0.
  236. assert_raises(ValueError, kmeans, X, 0)
  237. assert_raises(ValueError, kmeans2, X, 0)
  238. assert_raises(ValueError, kmeans2, X, np.array([]))
  239. def test_kmeans_large_thres(self):
  240. # Regression test for gh-1774
  241. x = np.array([1, 2, 3, 4, 10], dtype=float)
  242. res = kmeans(x, 1, thresh=1e16)
  243. assert_allclose(res[0], np.array([4.]))
  244. assert_allclose(res[1], 2.3999999999999999)
  245. def test_kmeans2_kpp_low_dim(self):
  246. # Regression test for gh-11462
  247. prev_res = np.array([[-1.95266667, 0.898],
  248. [-3.153375, 3.3945]])
  249. np.random.seed(42)
  250. res, _ = kmeans2(TESTDATA_2D, 2, minit='++')
  251. assert_allclose(res, prev_res)
  252. def test_kmeans2_kpp_high_dim(self):
  253. # Regression test for gh-11462
  254. n_dim = 100
  255. size = 10
  256. centers = np.vstack([5 * np.ones(n_dim),
  257. -5 * np.ones(n_dim)])
  258. np.random.seed(42)
  259. data = np.vstack([
  260. np.random.multivariate_normal(centers[0], np.eye(n_dim), size=size),
  261. np.random.multivariate_normal(centers[1], np.eye(n_dim), size=size)
  262. ])
  263. res, _ = kmeans2(data, 2, minit='++')
  264. assert_array_almost_equal(res, centers, decimal=0)
  265. def test_kmeans_and_kmeans2_random_seed(self):
  266. seed_list = [1234, np.random.RandomState(1234)]
  267. # check that np.random.Generator can be used (numpy >= 1.17)
  268. if hasattr(np.random, 'default_rng'):
  269. seed_list.append(np.random.default_rng(1234))
  270. for seed in seed_list:
  271. # test for kmeans
  272. res1, _ = kmeans(TESTDATA_2D, 2, seed=seed)
  273. res2, _ = kmeans(TESTDATA_2D, 2, seed=seed)
  274. assert_allclose(res1, res1) # should be same results
  275. # test for kmeans2
  276. for minit in ["random", "points", "++"]:
  277. res1, _ = kmeans2(TESTDATA_2D, 2, minit=minit, seed=seed)
  278. res2, _ = kmeans2(TESTDATA_2D, 2, minit=minit, seed=seed)
  279. assert_allclose(res1, res1) # should be same results