GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
find.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2015 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include "quit.h"
28 
29 #include "defun.h"
30 #include "error.h"
31 #include "gripes.h"
32 #include "oct-obj.h"
33 
34 // Find at most N_TO_FIND nonzero elements in NDA. Search forward if
35 // DIRECTION is 1, backward if it is -1. NARGOUT is the number of
36 // output arguments. If N_TO_FIND is -1, find all nonzero elements.
37 
38 template <typename T>
40 find_nonzero_elem_idx (const Array<T>& nda, int nargout,
41  octave_idx_type n_to_find, int direction)
42 {
43  octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ());
44 
46  if (n_to_find >= 0)
47  idx = nda.find (n_to_find, direction == -1);
48  else
49  idx = nda.find ();
50 
51  // The maximum element is always at the end.
52  octave_idx_type iext = idx.is_empty () ? 0 : idx.xelem (idx.numel () - 1) + 1;
53 
54  switch (nargout)
55  {
56  default:
57  case 3:
58  retval(2) = Array<T> (nda.index (idx_vector (idx)));
59  // Fall through!
60 
61  case 2:
62  {
63  Array<octave_idx_type> jdx (idx.dims ());
64  octave_idx_type n = idx.length ();
65  octave_idx_type nr = nda.rows ();
66  for (octave_idx_type i = 0; i < n; i++)
67  {
68  jdx.xelem (i) = idx.xelem (i) / nr;
69  idx.xelem (i) %= nr;
70  }
71  iext = -1;
72  retval(1) = idx_vector (jdx, -1);
73  }
74  // Fall through!
75 
76  case 1:
77  case 0:
78  retval(0) = idx_vector (idx, iext);
79  break;
80  }
81 
82  return retval;
83 }
84 
85 template <typename T>
87 find_nonzero_elem_idx (const Sparse<T>& v, int nargout,
88  octave_idx_type n_to_find, int direction)
89 {
90  octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ());
91 
92  octave_idx_type nr = v.rows ();
93  octave_idx_type nc = v.cols ();
94  octave_idx_type nz = v.nnz ();
95 
96  // Search in the default range.
97  octave_idx_type start_nc = -1;
98  octave_idx_type end_nc = -1;
99  octave_idx_type count;
100 
101  // Search for the range to search
102  if (n_to_find < 0)
103  {
104  start_nc = 0;
105  end_nc = nc;
106  n_to_find = nz;
107  count = nz;
108  }
109  else if (direction > 0)
110  {
111  for (octave_idx_type j = 0; j < nc; j++)
112  {
113  OCTAVE_QUIT;
114  if (v.cidx (j) == 0 && v.cidx (j+1) != 0)
115  start_nc = j;
116  if (v.cidx (j+1) >= n_to_find)
117  {
118  end_nc = j + 1;
119  break;
120  }
121  }
122  }
123  else
124  {
125  for (octave_idx_type j = nc; j > 0; j--)
126  {
127  OCTAVE_QUIT;
128  if (v.cidx (j) == nz && v.cidx (j-1) != nz)
129  end_nc = j;
130  if (nz - v.cidx (j-1) >= n_to_find)
131  {
132  start_nc = j - 1;
133  break;
134  }
135  }
136  }
137 
138  count = (n_to_find > v.cidx (end_nc) - v.cidx (start_nc) ?
139  v.cidx (end_nc) - v.cidx (start_nc) : n_to_find);
140 
141  octave_idx_type result_nr;
142  octave_idx_type result_nc;
143 
144  // Default case is to return a column vector, however, if the original
145  // argument was a row vector, then force return of a row vector.
146  if (nr == 1)
147  {
148  result_nr = 1;
149  result_nc = count;
150  }
151  else
152  {
153  result_nr = count;
154  result_nc = 1;
155  }
156 
157  Matrix idx (result_nr, result_nc);
158 
159  Matrix i_idx (result_nr, result_nc);
160  Matrix j_idx (result_nr, result_nc);
161 
162  Array<T> val (dim_vector (result_nr, result_nc));
163 
164  if (count > 0)
165  {
166  // Search for elements to return. Only search the region where there
167  // are elements to be found using the count that we want to find.
168  for (octave_idx_type j = start_nc, cx = 0; j < end_nc; j++)
169  for (octave_idx_type i = v.cidx (j); i < v.cidx (j+1); i++)
170  {
171  OCTAVE_QUIT;
172  if (direction < 0 && i < nz - count)
173  continue;
174  i_idx(cx) = static_cast<double> (v.ridx (i) + 1);
175  j_idx(cx) = static_cast<double> (j + 1);
176  idx(cx) = j * nr + v.ridx (i) + 1;
177  val(cx) = v.data(i);
178  cx++;
179  if (cx == count)
180  break;
181  }
182  }
183  else
184  {
185  // No items found. Fixup return dimensions for Matlab compatibility.
186  // The behavior to match is documented in Array.cc (Array<T>::find).
187  if ((nr == 0 && nc == 0) || (nr == 1 && nc == 1))
188  {
189  idx.resize (0, 0);
190 
191  i_idx.resize (0, 0);
192  j_idx.resize (0, 0);
193 
194  val.resize (dim_vector (0, 0));
195  }
196  }
197 
198  switch (nargout)
199  {
200  case 0:
201  case 1:
202  retval(0) = idx;
203  break;
204 
205  case 5:
206  retval(4) = nc;
207  // Fall through
208 
209  case 4:
210  retval(3) = nr;
211  // Fall through
212 
213  case 3:
214  retval(2) = val;
215  // Fall through!
216 
217  case 2:
218  retval(1) = j_idx;
219  retval(0) = i_idx;
220  break;
221 
222  default:
223  panic_impossible ();
224  break;
225  }
226 
227  return retval;
228 }
229 
231 find_nonzero_elem_idx (const PermMatrix& v, int nargout,
232  octave_idx_type n_to_find, int direction)
233 {
234  // There are far fewer special cases to handle for a PermMatrix.
235  octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ());
236 
237  octave_idx_type nr = v.rows ();
238  octave_idx_type nc = v.cols ();
239  octave_idx_type start_nc, count;
240 
241  // Determine the range to search.
242  if (n_to_find < 0 || n_to_find >= nc)
243  {
244  start_nc = 0;
245  count = nc;
246  }
247  else if (direction > 0)
248  {
249  start_nc = 0;
250  count = n_to_find;
251  }
252  else
253  {
254  start_nc = nc - n_to_find;
255  count = n_to_find;
256  }
257 
258  Matrix idx (count, 1);
259  Matrix i_idx (count, 1);
260  Matrix j_idx (count, 1);
261  // Every value is 1.
262  Array<double> val (dim_vector (count, 1), 1.0);
263 
264  if (count > 0)
265  {
266  const Array<octave_idx_type>& p = v.col_perm_vec ();
267  for (octave_idx_type k = 0; k < count; k++)
268  {
269  OCTAVE_QUIT;
270  const octave_idx_type j = start_nc + k;
271  const octave_idx_type i = p(j);
272  i_idx(k) = static_cast<double> (1+i);
273  j_idx(k) = static_cast<double> (1+j);
274  idx(k) = j * nc + i + 1;
275  }
276  }
277  else
278  {
279  // FIXME: Is this case even possible? A scalar permutation matrix seems
280  // to devolve to a scalar full matrix, at least from the Octave command
281  // line. Perhaps this function could be called internally from C++ with
282  // such a matrix.
283  // No items found. Fixup return dimensions for Matlab compatibility.
284  // The behavior to match is documented in Array.cc (Array<T>::find).
285  if ((nr == 0 && nc == 0) || (nr == 1 && nc == 1))
286  {
287  idx.resize (0, 0);
288 
289  i_idx.resize (0, 0);
290  j_idx.resize (0, 0);
291 
292  val.resize (dim_vector (0, 0));
293  }
294  }
295 
296  switch (nargout)
297  {
298  case 0:
299  case 1:
300  retval(0) = idx;
301  break;
302 
303  case 5:
304  retval(4) = nc;
305  // Fall through
306 
307  case 4:
308  retval(3) = nc;
309  // Fall through
310 
311  case 3:
312  retval(2) = val;
313  // Fall through!
314 
315  case 2:
316  retval(1) = j_idx;
317  retval(0) = i_idx;
318  break;
319 
320  default:
321  panic_impossible ();
322  break;
323  }
324 
325  return retval;
326 }
327 
328 DEFUN (find, args, nargout,
329  "-*- texinfo -*-\n\
330 @deftypefn {Built-in Function} {@var{idx} =} find (@var{x})\n\
331 @deftypefnx {Built-in Function} {@var{idx} =} find (@var{x}, @var{n})\n\
332 @deftypefnx {Built-in Function} {@var{idx} =} find (@var{x}, @var{n}, @var{direction})\n\
333 @deftypefnx {Built-in Function} {[i, j] =} find (@dots{})\n\
334 @deftypefnx {Built-in Function} {[i, j, v] =} find (@dots{})\n\
335 Return a vector of indices of nonzero elements of a matrix, as a row if\n\
336 @var{x} is a row vector or as a column otherwise.\n\
337 \n\
338 To obtain a single index for each matrix element, Octave pretends that the\n\
339 columns of a matrix form one long vector (like Fortran arrays are stored).\n\
340 For example:\n\
341 \n\
342 @example\n\
343 @group\n\
344 find (eye (2))\n\
345  @result{} [ 1; 4 ]\n\
346 @end group\n\
347 @end example\n\
348 \n\
349 If two inputs are given, @var{n} indicates the maximum number of elements to\n\
350 find from the beginning of the matrix or vector.\n\
351 \n\
352 If three inputs are given, @var{direction} should be one of\n\
353 @qcode{\"first\"} or @qcode{\"last\"}, requesting only the first or last\n\
354 @var{n} indices, respectively. However, the indices are always returned in\n\
355 ascending order.\n\
356 \n\
357 If two outputs are requested, @code{find} returns the row and column\n\
358 indices of nonzero elements of a matrix. For example:\n\
359 \n\
360 @example\n\
361 @group\n\
362 [i, j] = find (2 * eye (2))\n\
363  @result{} i = [ 1; 2 ]\n\
364  @result{} j = [ 1; 2 ]\n\
365 @end group\n\
366 @end example\n\
367 \n\
368 If three outputs are requested, @code{find} also returns a vector\n\
369 containing the nonzero values. For example:\n\
370 \n\
371 @example\n\
372 @group\n\
373 [i, j, v] = find (3 * eye (2))\n\
374  @result{} i = [ 1; 2 ]\n\
375  @result{} j = [ 1; 2 ]\n\
376  @result{} v = [ 3; 3 ]\n\
377 @end group\n\
378 @end example\n\
379 \n\
380 Note that this function is particularly useful for sparse matrices, as\n\
381 it extracts the nonzero elements as vectors, which can then be used to\n\
382 create the original matrix. For example:\n\
383 \n\
384 @example\n\
385 @group\n\
386 sz = size (a);\n\
387 [i, j, v] = find (a);\n\
388 b = sparse (i, j, v, sz(1), sz(2));\n\
389 @end group\n\
390 @end example\n\
391 @seealso{nonzeros}\n\
392 @end deftypefn")
393 {
394  octave_value_list retval;
395 
396  int nargin = args.length ();
397 
398  if (nargin > 3 || nargin < 1)
399  {
400  print_usage ();
401  return retval;
402  }
403 
404  // Setup the default options.
405  octave_idx_type n_to_find = -1;
406  if (nargin > 1)
407  {
408  double val = args(1).scalar_value ();
409 
410  if (error_state || (val < 0 || (! xisinf (val) && val != xround (val))))
411  {
412  error ("find: N must be a non-negative integer");
413  return retval;
414  }
415  else if (! xisinf (val))
416  n_to_find = val;
417  }
418 
419  // Direction to do the searching (1 == forward, -1 == reverse).
420  int direction = 1;
421  if (nargin > 2)
422  {
423  direction = 0;
424 
425  std::string s_arg = args(2).string_value ();
426 
427  if (! error_state)
428  {
429  if (s_arg == "first")
430  direction = 1;
431  else if (s_arg == "last")
432  direction = -1;
433  }
434 
435  if (direction == 0)
436  {
437  error ("find: DIRECTION must be \"first\" or \"last\"");
438  return retval;
439  }
440  }
441 
442  octave_value arg = args(0);
443 
444  if (arg.is_bool_type ())
445  {
446  if (arg.is_sparse_type ())
447  {
449 
450  if (! error_state)
451  retval = find_nonzero_elem_idx (v, nargout,
452  n_to_find, direction);
453  }
454  else if (nargout <= 1 && n_to_find == -1 && direction == 1)
455  {
456  // This case is equivalent to extracting indices from a logical
457  // matrix. Try to reuse the possibly cached index vector.
458  retval(0) = arg.index_vector ().unmask ();
459  }
460  else
461  {
462  boolNDArray v = arg.bool_array_value ();
463 
464  if (! error_state)
465  retval = find_nonzero_elem_idx (v, nargout,
466  n_to_find, direction);
467  }
468  }
469  else if (arg.is_integer_type ())
470  {
471 #define DO_INT_BRANCH(INTT) \
472  else if (arg.is_ ## INTT ## _type ()) \
473  { \
474  INTT ## NDArray v = arg.INTT ## _array_value (); \
475  \
476  if (! error_state) \
477  retval = find_nonzero_elem_idx (v, nargout, \
478  n_to_find, direction);\
479  }
480 
481  if (false)
482  ;
483  DO_INT_BRANCH (int8)
484  DO_INT_BRANCH (int16)
485  DO_INT_BRANCH (int32)
486  DO_INT_BRANCH (int64)
487  DO_INT_BRANCH (uint8)
488  DO_INT_BRANCH (uint16)
489  DO_INT_BRANCH (uint32)
490  DO_INT_BRANCH (uint64)
491  else
492  panic_impossible ();
493  }
494  else if (arg.is_sparse_type ())
495  {
496  if (arg.is_real_type ())
497  {
499 
500  if (! error_state)
501  retval = find_nonzero_elem_idx (v, nargout,
502  n_to_find, direction);
503  }
504  else if (arg.is_complex_type ())
505  {
507 
508  if (! error_state)
509  retval = find_nonzero_elem_idx (v, nargout,
510  n_to_find, direction);
511  }
512  else
513  gripe_wrong_type_arg ("find", arg);
514  }
515  else if (arg.is_perm_matrix ())
516  {
517  PermMatrix P = arg.perm_matrix_value ();
518 
519  if (! error_state)
520  retval = find_nonzero_elem_idx (P, nargout, n_to_find, direction);
521  }
522  else if (arg.is_string ())
523  {
524  charNDArray chnda = arg.char_array_value ();
525 
526  if (! error_state)
527  retval = find_nonzero_elem_idx (chnda, nargout, n_to_find, direction);
528  }
529  else if (arg.is_single_type ())
530  {
531  if (arg.is_real_type ())
532  {
533  FloatNDArray nda = arg.float_array_value ();
534 
535  if (! error_state)
536  retval = find_nonzero_elem_idx (nda, nargout, n_to_find,
537  direction);
538  }
539  else if (arg.is_complex_type ())
540  {
542 
543  if (! error_state)
544  retval = find_nonzero_elem_idx (cnda, nargout, n_to_find,
545  direction);
546  }
547  }
548  else if (arg.is_real_type ())
549  {
550  NDArray nda = arg.array_value ();
551 
552  if (! error_state)
553  retval = find_nonzero_elem_idx (nda, nargout, n_to_find, direction);
554  }
555  else if (arg.is_complex_type ())
556  {
557  ComplexNDArray cnda = arg.complex_array_value ();
558 
559  if (! error_state)
560  retval = find_nonzero_elem_idx (cnda, nargout, n_to_find, direction);
561  }
562  else
563  gripe_wrong_type_arg ("find", arg);
564 
565  return retval;
566 }
567 
568 /*
569 %!assert (find (char ([0, 97])), 2)
570 %!assert (find ([1, 0, 1, 0, 1]), [1, 3, 5])
571 %!assert (find ([1; 0; 3; 0; 1]), [1; 3; 5])
572 %!assert (find ([0, 0, 2; 0, 3, 0; -1, 0, 0]), [3; 5; 7])
573 
574 %!test
575 %! [i, j, v] = find ([0, 0, 2; 0, 3, 0; -1, 0, 0]);
576 %!
577 %! assert (i, [3; 2; 1]);
578 %! assert (j, [1; 2; 3]);
579 %! assert (v, [-1; 3; 2]);
580 
581 %!assert (find (single ([1, 0, 1, 0, 1])), [1, 3, 5])
582 %!assert (find (single ([1; 0; 3; 0; 1])), [1; 3; 5])
583 %!assert (find (single ([0, 0, 2; 0, 3, 0; -1, 0, 0])), [3; 5; 7])
584 
585 %!test
586 %! [i, j, v] = find (single ([0, 0, 2; 0, 3, 0; -1, 0, 0]));
587 %!
588 %! assert (i, [3; 2; 1]);
589 %! assert (j, [1; 2; 3]);
590 %! assert (v, single ([-1; 3; 2]));
591 
592 %!test
593 %! pcol = [5 1 4 3 2];
594 %! P = eye (5) (:, pcol);
595 %! [i, j, v] = find (P);
596 %! [ifull, jfull, vfull] = find (full (P));
597 %! assert (i, ifull);
598 %! assert (j, jfull);
599 %! assert (all (v == 1));
600 
601 %!test
602 %! prow = [5 1 4 3 2];
603 %! P = eye (5) (prow, :);
604 %! [i, j, v] = find (P);
605 %! [ifull, jfull, vfull] = find (full (P));
606 %! assert (i, ifull);
607 %! assert (j, jfull);
608 %! assert (all (v == 1));
609 
610 %!assert (find ([2 0 1 0 5 0], 1), 1)
611 %!assert (find ([2 0 1 0 5 0], 2, "last"), [3, 5])
612 
613 %!assert (find ([2 0 1 0 5 0], Inf), [1, 3, 5])
614 %!assert (find ([2 0 1 0 5 0], Inf, "last"), [1, 3, 5])
615 
616 %!error find ()
617 */
ComplexNDArray complex_array_value(bool frc_str_conv=false) const
Definition: ov.h:798
bool is_empty(void) const
Definition: Array.h:472
octave_value_list find_nonzero_elem_idx(const Array< T > &nda, int nargout, octave_idx_type n_to_find, int direction)
Definition: find.cc:40
octave_idx_type cols(void) const
Definition: Sparse.h:264
charNDArray char_array_value(bool frc_str_conv=false) const
Definition: ov.h:817
bool is_real_type(void) const
Definition: ov.h:651
T * data(void)
Definition: Sparse.h:509
void gripe_wrong_type_arg(const char *name, const char *s, bool is_error)
Definition: gripes.cc:135
octave_idx_type rows(void) const
Definition: Sparse.h:263
double xround(double x)
Definition: lo-mappers.cc:63
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:130
octave_idx_type rows(void) const
Definition: PermMatrix.h:55
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:275
octave_idx_type length(void) const
Definition: oct-obj.h:89
octave_idx_type cols(void) const
Definition: PermMatrix.h:56
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
bool is_perm_matrix(void) const
Definition: ov.h:559
bool xisinf(double x)
Definition: lo-mappers.cc:160
octave_idx_type * cidx(void)
Definition: Sparse.h:531
boolNDArray bool_array_value(bool warn=false) const
Definition: ov.h:811
octave_idx_type rows(void) const
Definition: Array.h:313
octave_idx_type nnz(void) const
Definition: Sparse.h:248
idx_vector index_vector(bool require_integers=false) const
Definition: ov.h:463
const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
Definition: Array.h:337
FloatNDArray float_array_value(bool frc_str_conv=false) const
Definition: ov.h:782
idx_vector unmask(void) const
Definition: idx-vector.cc:1209
bool is_sparse_type(void) const
Definition: ov.h:666
bool is_bool_type(void) const
Definition: ov.h:645
const Array< octave_idx_type > & col_perm_vec(void) const
Definition: PermMatrix.h:72
FloatComplexNDArray float_complex_array_value(bool frc_str_conv=false) const
Definition: ov.h:802
bool is_string(void) const
Definition: ov.h:562
int error_state
Definition: error.cc:101
SparseBoolMatrix sparse_bool_matrix_value(bool warn=false) const
Definition: ov.h:827
bool is_complex_type(void) const
Definition: ov.h:654
void resize(const dim_vector &dv, const T &rfv)
Definition: Array.cc:1033
#define panic_impossible()
Definition: error.h:33
Definition: dMatrix.h:35
SparseComplexMatrix sparse_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:824
#define DO_INT_BRANCH(INTT)
double arg(double x)
Definition: lo-mappers.h:37
static octave_idx_type find(octave_idx_type i, octave_idx_type *pp)
Definition: colamd.cc:111
T & xelem(octave_idx_type n)
Definition: Array.h:353
Handles the reference counting for all the derived classes.
Definition: Array.h:45
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
NDArray array_value(bool frc_str_conv=false) const
Definition: ov.h:779
octave_idx_type * ridx(void)
Definition: Sparse.h:518
Array< octave_idx_type > find(octave_idx_type n=-1, bool backward=false) const
Find indices of (at most n) nonzero elements.
Definition: Array.cc:2246
PermMatrix perm_matrix_value(void) const
Definition: ov.h:843
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:820
bool is_single_type(void) const
Definition: ov.h:611
#define OCTAVE_QUIT
Definition: quit.h:130
Array< T > index(const idx_vector &i) const
Indexing without resizing.
Definition: Array.cc:716
bool is_integer_type(void) const
Definition: ov.h:648