cxx11_tensor_reduction.cpp 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532
  1. // This file is part of Eigen, a lightweight C++ template library
  2. // for linear algebra.
  3. //
  4. // Copyright (C) 2014 Benoit Steiner <benoit.steiner.goog@gmail.com>
  5. //
  6. // This Source Code Form is subject to the terms of the Mozilla
  7. // Public License v. 2.0. If a copy of the MPL was not distributed
  8. // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
  9. #include "main.h"
  10. #include <limits>
  11. #include <numeric>
  12. #include <Eigen/CXX11/Tensor>
  13. using Eigen::Tensor;
  14. template <int DataLayout>
  15. static void test_trivial_reductions() {
  16. {
  17. Tensor<float, 0, DataLayout> tensor;
  18. tensor.setRandom();
  19. array<ptrdiff_t, 0> reduction_axis;
  20. Tensor<float, 0, DataLayout> result = tensor.sum(reduction_axis);
  21. VERIFY_IS_EQUAL(result(), tensor());
  22. }
  23. {
  24. Tensor<float, 1, DataLayout> tensor(7);
  25. tensor.setRandom();
  26. array<ptrdiff_t, 0> reduction_axis;
  27. Tensor<float, 1, DataLayout> result = tensor.sum(reduction_axis);
  28. VERIFY_IS_EQUAL(result.dimension(0), 7);
  29. for (int i = 0; i < 7; ++i) {
  30. VERIFY_IS_EQUAL(result(i), tensor(i));
  31. }
  32. }
  33. {
  34. Tensor<float, 2, DataLayout> tensor(2, 3);
  35. tensor.setRandom();
  36. array<ptrdiff_t, 0> reduction_axis;
  37. Tensor<float, 2, DataLayout> result = tensor.sum(reduction_axis);
  38. VERIFY_IS_EQUAL(result.dimension(0), 2);
  39. VERIFY_IS_EQUAL(result.dimension(1), 3);
  40. for (int i = 0; i < 2; ++i) {
  41. for (int j = 0; j < 3; ++j) {
  42. VERIFY_IS_EQUAL(result(i, j), tensor(i, j));
  43. }
  44. }
  45. }
  46. }
  47. template <typename Scalar,int DataLayout>
  48. static void test_simple_reductions() {
  49. Tensor<Scalar, 4, DataLayout> tensor(2, 3, 5, 7);
  50. tensor.setRandom();
  51. // Add a little offset so that the product reductions won't be close to zero.
  52. tensor += tensor.constant(Scalar(0.5f));
  53. array<ptrdiff_t, 2> reduction_axis2;
  54. reduction_axis2[0] = 1;
  55. reduction_axis2[1] = 3;
  56. Tensor<Scalar, 2, DataLayout> result = tensor.sum(reduction_axis2);
  57. VERIFY_IS_EQUAL(result.dimension(0), 2);
  58. VERIFY_IS_EQUAL(result.dimension(1), 5);
  59. for (int i = 0; i < 2; ++i) {
  60. for (int j = 0; j < 5; ++j) {
  61. Scalar sum = Scalar(0.0f);
  62. for (int k = 0; k < 3; ++k) {
  63. for (int l = 0; l < 7; ++l) {
  64. sum += tensor(i, k, j, l);
  65. }
  66. }
  67. VERIFY_IS_APPROX(result(i, j), sum);
  68. }
  69. }
  70. {
  71. Tensor<Scalar, 0, DataLayout> sum1 = tensor.sum();
  72. VERIFY_IS_EQUAL(sum1.rank(), 0);
  73. array<ptrdiff_t, 4> reduction_axis4;
  74. reduction_axis4[0] = 0;
  75. reduction_axis4[1] = 1;
  76. reduction_axis4[2] = 2;
  77. reduction_axis4[3] = 3;
  78. Tensor<Scalar, 0, DataLayout> sum2 = tensor.sum(reduction_axis4);
  79. VERIFY_IS_EQUAL(sum2.rank(), 0);
  80. VERIFY_IS_APPROX(sum1(), sum2());
  81. }
  82. reduction_axis2[0] = 0;
  83. reduction_axis2[1] = 2;
  84. result = tensor.prod(reduction_axis2);
  85. VERIFY_IS_EQUAL(result.dimension(0), 3);
  86. VERIFY_IS_EQUAL(result.dimension(1), 7);
  87. for (int i = 0; i < 3; ++i) {
  88. for (int j = 0; j < 7; ++j) {
  89. Scalar prod = Scalar(1.0f);
  90. for (int k = 0; k < 2; ++k) {
  91. for (int l = 0; l < 5; ++l) {
  92. prod *= tensor(k, i, l, j);
  93. }
  94. }
  95. VERIFY_IS_APPROX(result(i, j), prod);
  96. }
  97. }
  98. {
  99. Tensor<Scalar, 0, DataLayout> prod1 = tensor.prod();
  100. VERIFY_IS_EQUAL(prod1.rank(), 0);
  101. array<ptrdiff_t, 4> reduction_axis4;
  102. reduction_axis4[0] = 0;
  103. reduction_axis4[1] = 1;
  104. reduction_axis4[2] = 2;
  105. reduction_axis4[3] = 3;
  106. Tensor<Scalar, 0, DataLayout> prod2 = tensor.prod(reduction_axis4);
  107. VERIFY_IS_EQUAL(prod2.rank(), 0);
  108. VERIFY_IS_APPROX(prod1(), prod2());
  109. }
  110. reduction_axis2[0] = 0;
  111. reduction_axis2[1] = 2;
  112. result = tensor.maximum(reduction_axis2);
  113. VERIFY_IS_EQUAL(result.dimension(0), 3);
  114. VERIFY_IS_EQUAL(result.dimension(1), 7);
  115. for (int i = 0; i < 3; ++i) {
  116. for (int j = 0; j < 7; ++j) {
  117. Scalar max_val = std::numeric_limits<Scalar>::lowest();
  118. for (int k = 0; k < 2; ++k) {
  119. for (int l = 0; l < 5; ++l) {
  120. max_val = (std::max)(max_val, tensor(k, i, l, j));
  121. }
  122. }
  123. VERIFY_IS_APPROX(result(i, j), max_val);
  124. }
  125. }
  126. {
  127. Tensor<Scalar, 0, DataLayout> max1 = tensor.maximum();
  128. VERIFY_IS_EQUAL(max1.rank(), 0);
  129. array<ptrdiff_t, 4> reduction_axis4;
  130. reduction_axis4[0] = 0;
  131. reduction_axis4[1] = 1;
  132. reduction_axis4[2] = 2;
  133. reduction_axis4[3] = 3;
  134. Tensor<Scalar, 0, DataLayout> max2 = tensor.maximum(reduction_axis4);
  135. VERIFY_IS_EQUAL(max2.rank(), 0);
  136. VERIFY_IS_APPROX(max1(), max2());
  137. }
  138. reduction_axis2[0] = 0;
  139. reduction_axis2[1] = 1;
  140. result = tensor.minimum(reduction_axis2);
  141. VERIFY_IS_EQUAL(result.dimension(0), 5);
  142. VERIFY_IS_EQUAL(result.dimension(1), 7);
  143. for (int i = 0; i < 5; ++i) {
  144. for (int j = 0; j < 7; ++j) {
  145. Scalar min_val = (std::numeric_limits<Scalar>::max)();
  146. for (int k = 0; k < 2; ++k) {
  147. for (int l = 0; l < 3; ++l) {
  148. min_val = (std::min)(min_val, tensor(k, l, i, j));
  149. }
  150. }
  151. VERIFY_IS_APPROX(result(i, j), min_val);
  152. }
  153. }
  154. {
  155. Tensor<Scalar, 0, DataLayout> min1 = tensor.minimum();
  156. VERIFY_IS_EQUAL(min1.rank(), 0);
  157. array<ptrdiff_t, 4> reduction_axis4;
  158. reduction_axis4[0] = 0;
  159. reduction_axis4[1] = 1;
  160. reduction_axis4[2] = 2;
  161. reduction_axis4[3] = 3;
  162. Tensor<Scalar, 0, DataLayout> min2 = tensor.minimum(reduction_axis4);
  163. VERIFY_IS_EQUAL(min2.rank(), 0);
  164. VERIFY_IS_APPROX(min1(), min2());
  165. }
  166. reduction_axis2[0] = 0;
  167. reduction_axis2[1] = 1;
  168. result = tensor.mean(reduction_axis2);
  169. VERIFY_IS_EQUAL(result.dimension(0), 5);
  170. VERIFY_IS_EQUAL(result.dimension(1), 7);
  171. for (int i = 0; i < 5; ++i) {
  172. for (int j = 0; j < 7; ++j) {
  173. Scalar sum = Scalar(0.0f);
  174. int count = 0;
  175. for (int k = 0; k < 2; ++k) {
  176. for (int l = 0; l < 3; ++l) {
  177. sum += tensor(k, l, i, j);
  178. ++count;
  179. }
  180. }
  181. VERIFY_IS_APPROX(result(i, j), sum / Scalar(count));
  182. }
  183. }
  184. {
  185. Tensor<Scalar, 0, DataLayout> mean1 = tensor.mean();
  186. VERIFY_IS_EQUAL(mean1.rank(), 0);
  187. array<ptrdiff_t, 4> reduction_axis4;
  188. reduction_axis4[0] = 0;
  189. reduction_axis4[1] = 1;
  190. reduction_axis4[2] = 2;
  191. reduction_axis4[3] = 3;
  192. Tensor<Scalar, 0, DataLayout> mean2 = tensor.mean(reduction_axis4);
  193. VERIFY_IS_EQUAL(mean2.rank(), 0);
  194. VERIFY_IS_APPROX(mean1(), mean2());
  195. }
  196. {
  197. Tensor<int, 1> ints(10);
  198. std::iota(ints.data(), ints.data() + ints.dimension(0), 0);
  199. TensorFixedSize<bool, Sizes<> > all_;
  200. all_ = ints.all();
  201. VERIFY(!all_());
  202. all_ = (ints >= ints.constant(0)).all();
  203. VERIFY(all_());
  204. TensorFixedSize<bool, Sizes<> > any;
  205. any = (ints > ints.constant(10)).any();
  206. VERIFY(!any());
  207. any = (ints < ints.constant(1)).any();
  208. VERIFY(any());
  209. }
  210. }
  211. template <int DataLayout>
  212. static void test_reductions_in_expr() {
  213. Tensor<float, 4, DataLayout> tensor(2, 3, 5, 7);
  214. tensor.setRandom();
  215. array<ptrdiff_t, 2> reduction_axis2;
  216. reduction_axis2[0] = 1;
  217. reduction_axis2[1] = 3;
  218. Tensor<float, 2, DataLayout> result(2, 5);
  219. result = result.constant(1.0f) - tensor.sum(reduction_axis2);
  220. VERIFY_IS_EQUAL(result.dimension(0), 2);
  221. VERIFY_IS_EQUAL(result.dimension(1), 5);
  222. for (int i = 0; i < 2; ++i) {
  223. for (int j = 0; j < 5; ++j) {
  224. float sum = 0.0f;
  225. for (int k = 0; k < 3; ++k) {
  226. for (int l = 0; l < 7; ++l) {
  227. sum += tensor(i, k, j, l);
  228. }
  229. }
  230. VERIFY_IS_APPROX(result(i, j), 1.0f - sum);
  231. }
  232. }
  233. }
  234. template <int DataLayout>
  235. static void test_full_reductions() {
  236. Tensor<float, 2, DataLayout> tensor(2, 3);
  237. tensor.setRandom();
  238. array<ptrdiff_t, 2> reduction_axis;
  239. reduction_axis[0] = 0;
  240. reduction_axis[1] = 1;
  241. Tensor<float, 0, DataLayout> result = tensor.sum(reduction_axis);
  242. VERIFY_IS_EQUAL(result.rank(), 0);
  243. float sum = 0.0f;
  244. for (int i = 0; i < 2; ++i) {
  245. for (int j = 0; j < 3; ++j) {
  246. sum += tensor(i, j);
  247. }
  248. }
  249. VERIFY_IS_APPROX(result(0), sum);
  250. result = tensor.square().sum(reduction_axis).sqrt();
  251. VERIFY_IS_EQUAL(result.rank(), 0);
  252. sum = 0.0f;
  253. for (int i = 0; i < 2; ++i) {
  254. for (int j = 0; j < 3; ++j) {
  255. sum += tensor(i, j) * tensor(i, j);
  256. }
  257. }
  258. VERIFY_IS_APPROX(result(), sqrtf(sum));
  259. }
  260. struct UserReducer {
  261. static const bool PacketAccess = false;
  262. UserReducer(float offset) : offset_(offset) {}
  263. void reduce(const float val, float* accum) { *accum += val * val; }
  264. float initialize() const { return 0; }
  265. float finalize(const float accum) const { return 1.0f / (accum + offset_); }
  266. private:
  267. const float offset_;
  268. };
  269. template <int DataLayout>
  270. static void test_user_defined_reductions() {
  271. Tensor<float, 2, DataLayout> tensor(5, 7);
  272. tensor.setRandom();
  273. array<ptrdiff_t, 1> reduction_axis;
  274. reduction_axis[0] = 1;
  275. UserReducer reducer(10.0f);
  276. Tensor<float, 1, DataLayout> result = tensor.reduce(reduction_axis, reducer);
  277. VERIFY_IS_EQUAL(result.dimension(0), 5);
  278. for (int i = 0; i < 5; ++i) {
  279. float expected = 10.0f;
  280. for (int j = 0; j < 7; ++j) {
  281. expected += tensor(i, j) * tensor(i, j);
  282. }
  283. expected = 1.0f / expected;
  284. VERIFY_IS_APPROX(result(i), expected);
  285. }
  286. }
  287. template <int DataLayout>
  288. static void test_tensor_maps() {
  289. int inputs[2 * 3 * 5 * 7];
  290. TensorMap<Tensor<int, 4, DataLayout> > tensor_map(inputs, 2, 3, 5, 7);
  291. TensorMap<Tensor<const int, 4, DataLayout> > tensor_map_const(inputs, 2, 3, 5,
  292. 7);
  293. const TensorMap<Tensor<const int, 4, DataLayout> > tensor_map_const_const(
  294. inputs, 2, 3, 5, 7);
  295. tensor_map.setRandom();
  296. array<ptrdiff_t, 2> reduction_axis;
  297. reduction_axis[0] = 1;
  298. reduction_axis[1] = 3;
  299. Tensor<int, 2, DataLayout> result = tensor_map.sum(reduction_axis);
  300. Tensor<int, 2, DataLayout> result2 = tensor_map_const.sum(reduction_axis);
  301. Tensor<int, 2, DataLayout> result3 =
  302. tensor_map_const_const.sum(reduction_axis);
  303. for (int i = 0; i < 2; ++i) {
  304. for (int j = 0; j < 5; ++j) {
  305. int sum = 0;
  306. for (int k = 0; k < 3; ++k) {
  307. for (int l = 0; l < 7; ++l) {
  308. sum += tensor_map(i, k, j, l);
  309. }
  310. }
  311. VERIFY_IS_EQUAL(result(i, j), sum);
  312. VERIFY_IS_EQUAL(result2(i, j), sum);
  313. VERIFY_IS_EQUAL(result3(i, j), sum);
  314. }
  315. }
  316. }
  317. template <int DataLayout>
  318. static void test_static_dims() {
  319. Tensor<float, 4, DataLayout> in(72, 53, 97, 113);
  320. Tensor<float, 2, DataLayout> out(72, 97);
  321. in.setRandom();
  322. #if !EIGEN_HAS_CONSTEXPR
  323. array<int, 2> reduction_axis;
  324. reduction_axis[0] = 1;
  325. reduction_axis[1] = 3;
  326. #else
  327. Eigen::IndexList<Eigen::type2index<1>, Eigen::type2index<3> > reduction_axis;
  328. #endif
  329. out = in.maximum(reduction_axis);
  330. for (int i = 0; i < 72; ++i) {
  331. for (int j = 0; j < 97; ++j) {
  332. float expected = -1e10f;
  333. for (int k = 0; k < 53; ++k) {
  334. for (int l = 0; l < 113; ++l) {
  335. expected = (std::max)(expected, in(i, k, j, l));
  336. }
  337. }
  338. VERIFY_IS_EQUAL(out(i, j), expected);
  339. }
  340. }
  341. }
  342. template <int DataLayout>
  343. static void test_innermost_last_dims() {
  344. Tensor<float, 4, DataLayout> in(72, 53, 97, 113);
  345. Tensor<float, 2, DataLayout> out(97, 113);
  346. in.setRandom();
  347. // Reduce on the innermost dimensions.
  348. #if !EIGEN_HAS_CONSTEXPR
  349. array<int, 2> reduction_axis;
  350. reduction_axis[0] = 0;
  351. reduction_axis[1] = 1;
  352. #else
  353. // This triggers the use of packets for ColMajor.
  354. Eigen::IndexList<Eigen::type2index<0>, Eigen::type2index<1> > reduction_axis;
  355. #endif
  356. out = in.maximum(reduction_axis);
  357. for (int i = 0; i < 97; ++i) {
  358. for (int j = 0; j < 113; ++j) {
  359. float expected = -1e10f;
  360. for (int k = 0; k < 53; ++k) {
  361. for (int l = 0; l < 72; ++l) {
  362. expected = (std::max)(expected, in(l, k, i, j));
  363. }
  364. }
  365. VERIFY_IS_EQUAL(out(i, j), expected);
  366. }
  367. }
  368. }
  369. template <int DataLayout>
  370. static void test_innermost_first_dims() {
  371. Tensor<float, 4, DataLayout> in(72, 53, 97, 113);
  372. Tensor<float, 2, DataLayout> out(72, 53);
  373. in.setRandom();
  374. // Reduce on the innermost dimensions.
  375. #if !EIGEN_HAS_CONSTEXPR
  376. array<int, 2> reduction_axis;
  377. reduction_axis[0] = 2;
  378. reduction_axis[1] = 3;
  379. #else
  380. // This triggers the use of packets for RowMajor.
  381. Eigen::IndexList<Eigen::type2index<2>, Eigen::type2index<3>> reduction_axis;
  382. #endif
  383. out = in.maximum(reduction_axis);
  384. for (int i = 0; i < 72; ++i) {
  385. for (int j = 0; j < 53; ++j) {
  386. float expected = -1e10f;
  387. for (int k = 0; k < 97; ++k) {
  388. for (int l = 0; l < 113; ++l) {
  389. expected = (std::max)(expected, in(i, j, k, l));
  390. }
  391. }
  392. VERIFY_IS_EQUAL(out(i, j), expected);
  393. }
  394. }
  395. }
  396. template <int DataLayout>
  397. static void test_reduce_middle_dims() {
  398. Tensor<float, 4, DataLayout> in(72, 53, 97, 113);
  399. Tensor<float, 2, DataLayout> out(72, 53);
  400. in.setRandom();
  401. // Reduce on the innermost dimensions.
  402. #if !EIGEN_HAS_CONSTEXPR
  403. array<int, 2> reduction_axis;
  404. reduction_axis[0] = 1;
  405. reduction_axis[1] = 2;
  406. #else
  407. // This triggers the use of packets for RowMajor.
  408. Eigen::IndexList<Eigen::type2index<1>, Eigen::type2index<2>> reduction_axis;
  409. #endif
  410. out = in.maximum(reduction_axis);
  411. for (int i = 0; i < 72; ++i) {
  412. for (int j = 0; j < 113; ++j) {
  413. float expected = -1e10f;
  414. for (int k = 0; k < 53; ++k) {
  415. for (int l = 0; l < 97; ++l) {
  416. expected = (std::max)(expected, in(i, k, l, j));
  417. }
  418. }
  419. VERIFY_IS_EQUAL(out(i, j), expected);
  420. }
  421. }
  422. }
  423. static void test_sum_accuracy() {
  424. Tensor<float, 3> tensor(101, 101, 101);
  425. for (float prescribed_mean : {1.0f, 10.0f, 100.0f, 1000.0f, 10000.0f}) {
  426. tensor.setRandom();
  427. tensor += tensor.constant(prescribed_mean);
  428. Tensor<float, 0> sum = tensor.sum();
  429. double expected_sum = 0.0;
  430. for (int i = 0; i < 101; ++i) {
  431. for (int j = 0; j < 101; ++j) {
  432. for (int k = 0; k < 101; ++k) {
  433. expected_sum += static_cast<double>(tensor(i, j, k));
  434. }
  435. }
  436. }
  437. VERIFY_IS_APPROX(sum(), static_cast<float>(expected_sum));
  438. }
  439. }
  440. EIGEN_DECLARE_TEST(cxx11_tensor_reduction) {
  441. CALL_SUBTEST(test_trivial_reductions<ColMajor>());
  442. CALL_SUBTEST(test_trivial_reductions<RowMajor>());
  443. CALL_SUBTEST(( test_simple_reductions<float,ColMajor>() ));
  444. CALL_SUBTEST(( test_simple_reductions<float,RowMajor>() ));
  445. CALL_SUBTEST(( test_simple_reductions<Eigen::half,ColMajor>() ));
  446. CALL_SUBTEST(( test_simple_reductions<Eigen::bfloat16,ColMajor>() ));
  447. CALL_SUBTEST(test_reductions_in_expr<ColMajor>());
  448. CALL_SUBTEST(test_reductions_in_expr<RowMajor>());
  449. CALL_SUBTEST(test_full_reductions<ColMajor>());
  450. CALL_SUBTEST(test_full_reductions<RowMajor>());
  451. CALL_SUBTEST(test_user_defined_reductions<ColMajor>());
  452. CALL_SUBTEST(test_user_defined_reductions<RowMajor>());
  453. CALL_SUBTEST(test_tensor_maps<ColMajor>());
  454. CALL_SUBTEST(test_tensor_maps<RowMajor>());
  455. CALL_SUBTEST(test_static_dims<ColMajor>());
  456. CALL_SUBTEST(test_static_dims<RowMajor>());
  457. CALL_SUBTEST(test_innermost_last_dims<ColMajor>());
  458. CALL_SUBTEST(test_innermost_last_dims<RowMajor>());
  459. CALL_SUBTEST(test_innermost_first_dims<ColMajor>());
  460. CALL_SUBTEST(test_innermost_first_dims<RowMajor>());
  461. CALL_SUBTEST(test_reduce_middle_dims<ColMajor>());
  462. CALL_SUBTEST(test_reduce_middle_dims<RowMajor>());
  463. CALL_SUBTEST(test_sum_accuracy());
  464. }