spop_min_meat.hpp 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720
  1. // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
  2. // Copyright 2008-2016 National ICT Australia (NICTA)
  3. //
  4. // Licensed under the Apache License, Version 2.0 (the "License");
  5. // you may not use this file except in compliance with the License.
  6. // You may obtain a copy of the License at
  7. // http://www.apache.org/licenses/LICENSE-2.0
  8. //
  9. // Unless required by applicable law or agreed to in writing, software
  10. // distributed under the License is distributed on an "AS IS" BASIS,
  11. // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  12. // See the License for the specific language governing permissions and
  13. // limitations under the License.
  14. // ------------------------------------------------------------------------
  15. //! \addtogroup spop_min
  16. //! @{
  17. template<typename T1>
  18. inline
  19. void
  20. spop_min::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1,spop_min>& in)
  21. {
  22. arma_extra_debug_sigprint();
  23. const uword dim = in.aux_uword_a;
  24. arma_debug_check( (dim > 1), "min(): parameter 'dim' must be 0 or 1" );
  25. const SpProxy<T1> p(in.m);
  26. const uword p_n_rows = p.get_n_rows();
  27. const uword p_n_cols = p.get_n_cols();
  28. if( (p_n_rows == 0) || (p_n_cols == 0) || (p.get_n_nonzero() == 0) )
  29. {
  30. if(dim == 0) { out.zeros((p_n_rows > 0) ? 1 : 0, p_n_cols); }
  31. if(dim == 1) { out.zeros(p_n_rows, (p_n_cols > 0) ? 1 : 0); }
  32. return;
  33. }
  34. spop_min::apply_proxy(out, p, dim);
  35. }
  36. template<typename T1>
  37. inline
  38. void
  39. spop_min::apply_proxy
  40. (
  41. SpMat<typename T1::elem_type>& out,
  42. const SpProxy<T1>& p,
  43. const uword dim,
  44. const typename arma_not_cx<typename T1::elem_type>::result* junk
  45. )
  46. {
  47. arma_extra_debug_sigprint();
  48. arma_ignore(junk);
  49. typedef typename T1::elem_type eT;
  50. typename SpProxy<T1>::const_iterator_type it = p.begin();
  51. typename SpProxy<T1>::const_iterator_type it_end = p.end();
  52. const uword p_n_cols = p.get_n_cols();
  53. const uword p_n_rows = p.get_n_rows();
  54. if(dim == 0) // find the minimum in each column
  55. {
  56. Row<eT> value(p_n_cols, fill::zeros);
  57. urowvec count(p_n_cols, fill::zeros);
  58. while(it != it_end)
  59. {
  60. const uword col = it.col();
  61. value[col] = (count[col] == 0) ? (*it) : (std::min)(value[col], (*it));
  62. count[col]++;
  63. ++it;
  64. }
  65. for(uword col=0; col<p_n_cols; ++col)
  66. {
  67. if(count[col] < p_n_rows) { value[col] = (std::min)(value[col], eT(0)); }
  68. }
  69. out = value;
  70. }
  71. else
  72. if(dim == 1) // find the minimum in each row
  73. {
  74. Col<eT> value(p_n_rows, fill::zeros);
  75. ucolvec count(p_n_rows, fill::zeros);
  76. while(it != it_end)
  77. {
  78. const uword row = it.row();
  79. value[row] = (count[row] == 0) ? (*it) : (std::min)(value[row], (*it));
  80. count[row]++;
  81. ++it;
  82. }
  83. for(uword row=0; row<p_n_rows; ++row)
  84. {
  85. if(count[row] < p_n_cols) { value[row] = (std::min)(value[row], eT(0)); }
  86. }
  87. out = value;
  88. }
  89. }
  90. template<typename T1>
  91. inline
  92. typename T1::elem_type
  93. spop_min::vector_min
  94. (
  95. const T1& x,
  96. const typename arma_not_cx<typename T1::elem_type>::result* junk
  97. )
  98. {
  99. arma_extra_debug_sigprint();
  100. arma_ignore(junk);
  101. typedef typename T1::elem_type eT;
  102. const SpProxy<T1> p(x);
  103. if(p.get_n_elem() == 0)
  104. {
  105. arma_debug_check(true, "min(): object has no elements");
  106. return Datum<eT>::nan;
  107. }
  108. if(p.get_n_nonzero() == 0) { return eT(0); }
  109. if(SpProxy<T1>::use_iterator == false)
  110. {
  111. // direct access of values
  112. if(p.get_n_nonzero() == p.get_n_elem())
  113. {
  114. return op_min::direct_min(p.get_values(), p.get_n_nonzero());
  115. }
  116. else
  117. {
  118. return std::min(eT(0), op_min::direct_min(p.get_values(), p.get_n_nonzero()));
  119. }
  120. }
  121. else
  122. {
  123. // use iterator
  124. typename SpProxy<T1>::const_iterator_type it = p.begin();
  125. typename SpProxy<T1>::const_iterator_type it_end = p.end();
  126. eT result = (*it);
  127. ++it;
  128. while(it != it_end)
  129. {
  130. if((*it) < result) { result = (*it); }
  131. ++it;
  132. }
  133. if(p.get_n_nonzero() == p.get_n_elem())
  134. {
  135. return result;
  136. }
  137. else
  138. {
  139. return std::min(eT(0), result);
  140. }
  141. }
  142. }
  143. template<typename T1>
  144. inline
  145. typename arma_not_cx<typename T1::elem_type>::result
  146. spop_min::min(const SpBase<typename T1::elem_type, T1>& X)
  147. {
  148. arma_extra_debug_sigprint();
  149. typedef typename T1::elem_type eT;
  150. const SpProxy<T1> P(X.get_ref());
  151. const uword n_elem = P.get_n_elem();
  152. const uword n_nonzero = P.get_n_nonzero();
  153. if(n_elem == 0)
  154. {
  155. arma_debug_check(true, "min(): object has no elements");
  156. return Datum<eT>::nan;
  157. }
  158. eT min_val = priv::most_pos<eT>();
  159. if(SpProxy<T1>::use_iterator)
  160. {
  161. // We have to iterate over the elements.
  162. typedef typename SpProxy<T1>::const_iterator_type it_type;
  163. it_type it = P.begin();
  164. it_type it_end = P.end();
  165. while (it != it_end)
  166. {
  167. if ((*it) < min_val) { min_val = *it; }
  168. ++it;
  169. }
  170. }
  171. else
  172. {
  173. // We can do direct access of the values, row_indices, and col_ptrs.
  174. // We don't need the location of the min value, so we can just call out to
  175. // other functions...
  176. min_val = op_min::direct_min(P.get_values(), n_nonzero);
  177. }
  178. if(n_elem == n_nonzero)
  179. {
  180. return min_val;
  181. }
  182. else
  183. {
  184. return std::min(eT(0), min_val);
  185. }
  186. }
  187. template<typename T1>
  188. inline
  189. typename arma_not_cx<typename T1::elem_type>::result
  190. spop_min::min_with_index(const SpProxy<T1>& P, uword& index_of_min_val)
  191. {
  192. arma_extra_debug_sigprint();
  193. typedef typename T1::elem_type eT;
  194. const uword n_elem = P.get_n_elem();
  195. const uword n_nonzero = P.get_n_nonzero();
  196. const uword n_rows = P.get_n_rows();
  197. if(n_elem == 0)
  198. {
  199. arma_debug_check(true, "min(): object has no elements");
  200. index_of_min_val = uword(0);
  201. return Datum<eT>::nan;
  202. }
  203. eT min_val = priv::most_pos<eT>();
  204. if(SpProxy<T1>::use_iterator)
  205. {
  206. // We have to iterate over the elements.
  207. typedef typename SpProxy<T1>::const_iterator_type it_type;
  208. it_type it = P.begin();
  209. it_type it_end = P.end();
  210. while (it != it_end)
  211. {
  212. if ((*it) < min_val)
  213. {
  214. min_val = *it;
  215. index_of_min_val = it.row() + it.col() * n_rows;
  216. }
  217. ++it;
  218. }
  219. }
  220. else
  221. {
  222. // We can do direct access.
  223. min_val = op_min::direct_min(P.get_values(), n_nonzero, index_of_min_val);
  224. // Convert to actual position in matrix.
  225. const uword row = P.get_row_indices()[index_of_min_val];
  226. uword col = 0;
  227. while (P.get_col_ptrs()[++col] < index_of_min_val + 1) { }
  228. index_of_min_val = (col - 1) * n_rows + row;
  229. }
  230. if(n_elem != n_nonzero)
  231. {
  232. min_val = std::min(eT(0), min_val);
  233. // If the min_val is a nonzero element, we need its actual position in the matrix.
  234. if(min_val == eT(0))
  235. {
  236. // Find first zero element.
  237. uword last_row = 0;
  238. uword last_col = 0;
  239. typedef typename SpProxy<T1>::const_iterator_type it_type;
  240. it_type it = P.begin();
  241. it_type it_end = P.end();
  242. while (it != it_end)
  243. {
  244. // Have we moved more than one position from the last place?
  245. if ((it.col() == last_col) && (it.row() - last_row > 1))
  246. {
  247. index_of_min_val = it.col() * n_rows + last_row + 1;
  248. break;
  249. }
  250. else if ((it.col() >= last_col + 1) && (last_row < n_rows - 1))
  251. {
  252. index_of_min_val = last_col * n_rows + last_row + 1;
  253. break;
  254. }
  255. else if ((it.col() == last_col + 1) && (it.row() > 0))
  256. {
  257. index_of_min_val = it.col() * n_rows;
  258. break;
  259. }
  260. else if (it.col() > last_col + 1)
  261. {
  262. index_of_min_val = (last_col + 1) * n_rows;
  263. break;
  264. }
  265. last_row = it.row();
  266. last_col = it.col();
  267. ++it;
  268. }
  269. }
  270. }
  271. return min_val;
  272. }
  273. template<typename T1>
  274. inline
  275. void
  276. spop_min::apply_proxy
  277. (
  278. SpMat<typename T1::elem_type>& out,
  279. const SpProxy<T1>& p,
  280. const uword dim,
  281. const typename arma_cx_only<typename T1::elem_type>::result* junk
  282. )
  283. {
  284. arma_extra_debug_sigprint();
  285. arma_ignore(junk);
  286. typedef typename T1::elem_type eT;
  287. typedef typename get_pod_type<eT>::result T;
  288. typename SpProxy<T1>::const_iterator_type it = p.begin();
  289. typename SpProxy<T1>::const_iterator_type it_end = p.end();
  290. const uword p_n_cols = p.get_n_cols();
  291. const uword p_n_rows = p.get_n_rows();
  292. if(dim == 0) // find the minimum in each column
  293. {
  294. Row<eT> rawval(p_n_cols, fill::zeros);
  295. Row< T> absval(p_n_cols, fill::zeros);
  296. urowvec count(p_n_cols, fill::zeros);
  297. while(it != it_end)
  298. {
  299. const uword col = it.col();
  300. const eT& v = (*it);
  301. const T a = std::abs(v);
  302. if(count[col] == 0)
  303. {
  304. absval[col] = a;
  305. rawval[col] = v;
  306. }
  307. else
  308. {
  309. if(a < absval[col])
  310. {
  311. absval[col] = a;
  312. rawval[col] = v;
  313. }
  314. }
  315. count[col]++;
  316. ++it;
  317. }
  318. for(uword col=0; col < p_n_cols; ++col)
  319. {
  320. if(count[col] < p_n_rows)
  321. {
  322. if(T(0) < absval[col]) { rawval[col] = eT(0); }
  323. }
  324. }
  325. out = rawval;
  326. }
  327. else
  328. if(dim == 1) // find the minimum in each row
  329. {
  330. Col<eT> rawval(p_n_rows, fill::zeros);
  331. Col< T> absval(p_n_rows, fill::zeros);
  332. ucolvec count(p_n_rows, fill::zeros);
  333. while(it != it_end)
  334. {
  335. const uword row = it.row();
  336. const eT& v = (*it);
  337. const T a = std::abs(v);
  338. if(count[row] == 0)
  339. {
  340. absval[row] = a;
  341. rawval[row] = v;
  342. }
  343. else
  344. {
  345. if(a < absval[row])
  346. {
  347. absval[row] = a;
  348. rawval[row] = v;
  349. }
  350. }
  351. count[row]++;
  352. ++it;
  353. }
  354. for(uword row=0; row < p_n_rows; ++row)
  355. {
  356. if(count[row] < p_n_cols)
  357. {
  358. if(T(0) < absval[row]) { rawval[row] = eT(0); }
  359. }
  360. }
  361. out = rawval;
  362. }
  363. }
  364. template<typename T1>
  365. inline
  366. typename T1::elem_type
  367. spop_min::vector_min
  368. (
  369. const T1& x,
  370. const typename arma_cx_only<typename T1::elem_type>::result* junk
  371. )
  372. {
  373. arma_extra_debug_sigprint();
  374. arma_ignore(junk);
  375. typedef typename T1::elem_type eT;
  376. typedef typename get_pod_type<eT>::result T;
  377. const SpProxy<T1> p(x);
  378. if(p.get_n_elem() == 0)
  379. {
  380. arma_debug_check(true, "min(): object has no elements");
  381. return Datum<eT>::nan;
  382. }
  383. if(p.get_n_nonzero() == 0) { return eT(0); }
  384. if(SpProxy<T1>::use_iterator == false)
  385. {
  386. // direct access of values
  387. if(p.get_n_nonzero() == p.get_n_elem())
  388. {
  389. return op_min::direct_min(p.get_values(), p.get_n_nonzero());
  390. }
  391. else
  392. {
  393. const eT val1 = eT(0);
  394. const eT val2 = op_min::direct_min(p.get_values(), p.get_n_nonzero());
  395. return ( std::abs(val1) < std::abs(val2) ) ? val1 : val2;
  396. }
  397. }
  398. else
  399. {
  400. // use iterator
  401. typename SpProxy<T1>::const_iterator_type it = p.begin();
  402. typename SpProxy<T1>::const_iterator_type it_end = p.end();
  403. eT best_val_orig = *it;
  404. T best_val_abs = std::abs(best_val_orig);
  405. ++it;
  406. while(it != it_end)
  407. {
  408. eT val_orig = *it;
  409. T val_abs = std::abs(val_orig);
  410. if(val_abs < best_val_abs)
  411. {
  412. best_val_abs = val_abs;
  413. best_val_orig = val_orig;
  414. }
  415. ++it;
  416. }
  417. if(p.get_n_nonzero() == p.get_n_elem())
  418. {
  419. return best_val_orig;
  420. }
  421. else
  422. {
  423. const eT val1 = eT(0);
  424. return ( std::abs(val1) < best_val_abs ) ? val1 : best_val_orig;
  425. }
  426. }
  427. }
  428. template<typename T1>
  429. inline
  430. typename arma_cx_only<typename T1::elem_type>::result
  431. spop_min::min(const SpBase<typename T1::elem_type, T1>& X)
  432. {
  433. arma_extra_debug_sigprint();
  434. typedef typename T1::elem_type eT;
  435. typedef typename get_pod_type<eT>::result T;
  436. const SpProxy<T1> P(X.get_ref());
  437. const uword n_elem = P.get_n_elem();
  438. const uword n_nonzero = P.get_n_nonzero();
  439. if(n_elem == 0)
  440. {
  441. arma_debug_check(true, "min(): object has no elements");
  442. return Datum<eT>::nan;
  443. }
  444. T min_val = priv::most_pos<T>();
  445. eT ret_val;
  446. if(SpProxy<T1>::use_iterator)
  447. {
  448. // We have to iterate over the elements.
  449. typedef typename SpProxy<T1>::const_iterator_type it_type;
  450. it_type it = P.begin();
  451. it_type it_end = P.end();
  452. while (it != it_end)
  453. {
  454. const T tmp_val = std::abs(*it);
  455. if (tmp_val < min_val)
  456. {
  457. min_val = tmp_val;
  458. ret_val = *it;
  459. }
  460. ++it;
  461. }
  462. }
  463. else
  464. {
  465. // We can do direct access of the values, row_indices, and col_ptrs.
  466. // We don't need the location of the min value, so we can just call out to
  467. // other functions...
  468. ret_val = op_min::direct_min(P.get_values(), n_nonzero);
  469. min_val = std::abs(ret_val);
  470. }
  471. if(n_elem == n_nonzero)
  472. {
  473. return ret_val;
  474. }
  475. else
  476. {
  477. return (T(0) < min_val) ? eT(0) : ret_val;
  478. }
  479. }
  480. template<typename T1>
  481. inline
  482. typename arma_cx_only<typename T1::elem_type>::result
  483. spop_min::min_with_index(const SpProxy<T1>& P, uword& index_of_min_val)
  484. {
  485. arma_extra_debug_sigprint();
  486. typedef typename T1::elem_type eT;
  487. typedef typename get_pod_type<eT>::result T;
  488. const uword n_elem = P.get_n_elem();
  489. const uword n_nonzero = P.get_n_nonzero();
  490. const uword n_rows = P.get_n_rows();
  491. if(n_elem == 0)
  492. {
  493. arma_debug_check(true, "min(): object has no elements");
  494. index_of_min_val = uword(0);
  495. return Datum<eT>::nan;
  496. }
  497. T min_val = priv::most_pos<T>();
  498. if(SpProxy<T1>::use_iterator)
  499. {
  500. // We have to iterate over the elements.
  501. typedef typename SpProxy<T1>::const_iterator_type it_type;
  502. it_type it = P.begin();
  503. it_type it_end = P.end();
  504. while (it != it_end)
  505. {
  506. const T tmp_val = std::abs(*it);
  507. if (tmp_val < min_val)
  508. {
  509. min_val = tmp_val;
  510. index_of_min_val = it.row() + it.col() * n_rows;
  511. }
  512. ++it;
  513. }
  514. }
  515. else
  516. {
  517. // We can do direct access.
  518. min_val = std::abs(op_min::direct_min(P.get_values(), n_nonzero, index_of_min_val));
  519. // Convert to actual position in matrix.
  520. const uword row = P.get_row_indices()[index_of_min_val];
  521. uword col = 0;
  522. while (P.get_col_ptrs()[++col] < index_of_min_val + 1) { }
  523. index_of_min_val = (col - 1) * n_rows + row;
  524. }
  525. if(n_elem != n_nonzero)
  526. {
  527. min_val = std::min(T(0), min_val);
  528. // If the min_val is a nonzero element, we need its actual position in the matrix.
  529. if(min_val == T(0))
  530. {
  531. // Find first zero element.
  532. uword last_row = 0;
  533. uword last_col = 0;
  534. typedef typename SpProxy<T1>::const_iterator_type it_type;
  535. it_type it = P.begin();
  536. it_type it_end = P.end();
  537. while (it != it_end)
  538. {
  539. // Have we moved more than one position from the last place?
  540. if ((it.col() == last_col) && (it.row() - last_row > 1))
  541. {
  542. index_of_min_val = it.col() * n_rows + last_row + 1;
  543. break;
  544. }
  545. else if ((it.col() >= last_col + 1) && (last_row < n_rows - 1))
  546. {
  547. index_of_min_val = last_col * n_rows + last_row + 1;
  548. break;
  549. }
  550. else if ((it.col() == last_col + 1) && (it.row() > 0))
  551. {
  552. index_of_min_val = it.col() * n_rows;
  553. break;
  554. }
  555. else if (it.col() > last_col + 1)
  556. {
  557. index_of_min_val = (last_col + 1) * n_rows;
  558. break;
  559. }
  560. last_row = it.row();
  561. last_col = it.col();
  562. ++it;
  563. }
  564. }
  565. }
  566. return P[index_of_min_val];
  567. }
  568. //! @}