Mercurial > octave-nkf

changeset 9043:cfad8f9a77fa

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Documentation merge to main branch
author Rik <rdrider0-list@yahoo.com>
date Thu, 26 Mar 2009 09:47:35 -0700
parents 97aa01a85ea4 (current diff) 2df28ad88b0e (diff)
children 656ad518f385
files src/DLD-FUNCTIONS/find.cc
diffstat 9 files changed, 218 insertions(+), 142 deletions(-) [+]
line wrap: on
line diff
--- a/.hgtags	Wed Mar 25 18:13:08 2009 -0700
+++ b/.hgtags	Thu Mar 26 09:47:35 2009 -0700
@@ -41,3 +41,4 @@
 739141cde75a22b9ac9b5b327ea30c61fd357748 ss-3-1-52
 bd1b1fe9c6e952202bf5878cd042a168f6df6ddd ss-3-1-53
 3c3cbe8f18e03c73cc34494c8037d8a45a30bf7d ss-3-1-54
+5e276a0b999747b72af34e72eb63f64640da6399 ss-3-1-55
--- a/liboctave/Array.cc	Wed Mar 25 18:13:08 2009 -0700
+++ b/liboctave/Array.cc	Thu Mar 26 09:47:35 2009 -0700
@@ -2487,6 +2487,69 @@
   return idx;
 }
 
+template <class T>
+Array<octave_idx_type> 
+Array<T>::find (octave_idx_type n, bool backward) const
+{
+  Array<octave_idx_type> retval;
+  const T *src = data ();
+  octave_idx_type nel = nelem ();
+  const T zero = T ();
+  if (n < 0)
+    {
+      // We want all elements, which means we'll almost surely need
+      // to resize. So count first, then allocate array of exact size.
+      octave_idx_type cnt = 0;
+      for (octave_idx_type i = 0; i < nel; i++)
+        cnt += src[i] != zero;
+
+      retval = Array<octave_idx_type> (cnt);
+      octave_idx_type *dest = retval.fortran_vec ();
+      for (octave_idx_type i = 0; i < nel; i++)
+        if (src[i] != zero) *dest++ = i;
+    }
+  else
+    {
+      // We want a fixed max number of elements, usually small. So be
+      // optimistic, alloc the array in advance, and then resize if
+      // needed.
+      retval = Array<octave_idx_type> (n);
+      if (backward)
+        {
+          // Do the search as a series of successive single-element searches.
+          octave_idx_type k = 0, l = nel - 1;
+          for (; k < n; k++)
+            {
+              for (;l >= 0 && src[l] == zero; l--) ;
+              if (l >= 0)
+                retval(k) = l--;
+              else
+                break;
+            }
+          if (k < n)
+            retval.resize (k);
+        }
+      else
+        {
+          // Do the search as a series of successive single-element searches.
+          octave_idx_type k = 0, l = 0;
+          for (; k < n; k++)
+            {
+              for (;l != nel && src[l] == zero; l++) ;
+              if (l != nel)
+                retval(k) = l++;
+              else
+                break;
+            }
+          if (k < n)
+            retval.resize (k);
+        }
+    }
+
+  return retval;
+}
+
+
 #define INSTANTIATE_ARRAY_SORT(T) template class octave_sort<T>;
 
 #define NO_INSTANTIATE_ARRAY_SORT(T) \
@@ -2520,7 +2583,9 @@
 template <> Array<octave_idx_type>  \
 Array<T>::lookup (const Array<T>&, sortmode, bool, bool) const \
 { return Array<octave_idx_type> (); } \
-
+template <> Array<octave_idx_type> \
+Array<T>::find (octave_idx_type, bool) const\
+{ return Array<octave_idx_type> (); } \
 
 
 template <class T>
--- a/liboctave/Array.h	Wed Mar 25 18:13:08 2009 -0700
+++ b/liboctave/Array.h	Thu Mar 26 09:47:35 2009 -0700
@@ -255,7 +255,8 @@
 
   size_t byte_size (void) const { return numel () * sizeof (T); }
 
-  dim_vector dims (void) const { return dimensions; }
+  // Return a const-reference so that dims ()(i) works efficiently.
+  const dim_vector& dims (void) const { return dimensions; }
 
   Array<T> squeeze (void) const;
   
@@ -428,6 +429,8 @@
 
   bool is_empty (void) const { return numel () == 0; }
 
+  bool is_vector (void) const { return dimensions.is_vector (); }
+
   Array<T> transpose (void) const;
   Array<T> hermitian (T (*fcn) (const T&) = 0) const;
 
@@ -579,6 +582,10 @@
   Array<octave_idx_type> lookup (const Array<T>& values, sortmode mode = UNSORTED, 
                                  bool linf = false, bool rinf = false) const;
 
+  // Find indices of (at most n) nonzero elements. If n is specified, backward
+  // specifies search from backward.
+  Array<octave_idx_type> find (octave_idx_type n = -1, bool backward = false) const;
+
   Array<T> diag (octave_idx_type k = 0) const;
 
   template <class U, class F>
--- a/liboctave/ChangeLog	Wed Mar 25 18:13:08 2009 -0700
+++ b/liboctave/ChangeLog	Thu Mar 26 09:47:35 2009 -0700
@@ -1,5 +1,23 @@
+2009-03-26  Jaroslav Hajek  <highegg@gmail.com>
+
+	* dim-vector.h (dim_vector::numel): Add optional argument, simplify.
+
+2009-03-26  Jaroslav Hajek  <highegg@gmail.com>
+
+	* Array.h (Array<T>::dims): Return a const reference.
+	(Array<T>::is_vector): New method.
+
+2009-03-26  Jaroslav Hajek  <highegg@gmail.com>
+
+	* Array.cc (Array<T>::find): New method.
+	* Array.h: Declare it.
+
 2009-03-25  John W. Eaton  <jwe@octave.org>
 
+	* EIG.cc (EIG::init (const Matrix&, bool),
+	EIG::init (const Matrix&, const Matrix&, bool)):
+	Avoid volatile declaration for tmp variable.
+
 	* Makefile.in (MATRIX_INC): Add Sparse-diag-op-defs.h and
 	Sparse-perm-op-defs.h to the list.
 
--- a/liboctave/EIG.cc	Wed Mar 25 18:13:08 2009 -0700
+++ b/liboctave/EIG.cc	Thu Mar 26 09:47:35 2009 -0700
@@ -161,8 +161,8 @@
   Array<double> wi (n);
   double *pwi = wi.fortran_vec ();
 
-  volatile octave_idx_type nvr = calc_ev ? n : 0;
-  Matrix vr (nvr, nvr);
+  octave_idx_type tnvr = calc_ev ? n : 0;
+  Matrix vr (tnvr, tnvr);
   double *pvr = vr.fortran_vec ();
 
   octave_idx_type lwork = -1;
@@ -204,6 +204,7 @@
 	}
 
       lambda.resize (n);
+      octave_idx_type nvr = calc_ev ? n : 0;
       v.resize (nvr, nvr);
 
       for (octave_idx_type j = 0; j < n; j++)
@@ -507,8 +508,8 @@
   Array<double> beta (n);
   double *pbeta = beta.fortran_vec ();
 
-  volatile octave_idx_type nvr = calc_ev ? n : 0;
-  Matrix vr (nvr, nvr);
+  octave_idx_type tnvr = calc_ev ? n : 0;
+  Matrix vr (tnvr, tnvr);
   double *pvr = vr.fortran_vec ();
 
   octave_idx_type lwork = -1;
@@ -554,6 +555,7 @@
 	}
 
       lambda.resize (n);
+      octave_idx_type nvr = calc_ev ? n : 0;
       v.resize (nvr, nvr);
 
       for (octave_idx_type j = 0; j < n; j++)
--- a/liboctave/dim-vector.h	Wed Mar 25 18:13:08 2009 -0700
+++ b/liboctave/dim-vector.h	Thu Mar 26 09:47:35 2009 -0700
@@ -310,13 +310,13 @@
   // vector would have, NOT the number of dimensions (elements in the
   // dimension vector).
 
-  octave_idx_type numel (void) const
+  octave_idx_type numel (int n = 0) const
   {
     int n_dims = length ();
 
-    octave_idx_type retval = n_dims > 0 ? elem (0) : 0;
+    octave_idx_type retval = 1;
 
-    for (int i = 1; i < n_dims; i++)
+    for (int i = n; i < n_dims; i++)
       retval *= elem (i);
 
     return retval;
--- a/src/ChangeLog	Wed Mar 25 18:13:08 2009 -0700
+++ b/src/ChangeLog	Thu Mar 26 09:47:35 2009 -0700
@@ -1,5 +1,21 @@
+2009-03-26  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DLD-FUNCTIONS/find.cc 
+	(find_nonzero_elem_idx (const Array<T>&, ...)): Simplify.
+	Instantiate for bool and octave_int types.
+	(find_nonzero_elem_idx (const Sparse<T>&, ...)): 
+	Instantiate for bool.
+	(Ffind): Handle bool and octave_int cases.
+
 2009-03-25  John W. Eaton  <jwe@octave.org>
 
+	* version.h (OCTAVE_VERSION): Now 3.1.55+.
+	(OCTAVE_API_VERSION): Now api-v37+.
+
+	* version.h (OCTAVE_VERSION): Now 3.1.55.
+	(OCTAVE_API_VERSION): Now api-v37.
+	(OCTAVE_RELEASE_DATE): Now 2009-03-25.
+
 	* DLD-FUNCTIONS/find.cc (find_nonzero_elem_idx): Also return
 	[](0x0) if the array has 0 rows and it is not a column vector.
 
--- a/src/DLD-FUNCTIONS/find.cc	Wed Mar 25 18:13:08 2009 -0700
+++ b/src/DLD-FUNCTIONS/find.cc	Thu Mar 26 09:47:35 2009 -0700
@@ -43,155 +43,75 @@
 {
   octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ());
 
-  octave_idx_type count = 0;
-
-  octave_idx_type nel = nda.nelem ();
-
-  // Set the starting element to the correct value based on the
-  // direction to search.
-  octave_idx_type k = 0;
-  if (direction == -1)
-    k = nel - 1;
-
-  // Search in the default range.
-  octave_idx_type start_el = -1;
-  octave_idx_type end_el = -1;
-
-  // Search for the number of elements to return.
-  while (k < nel && k > -1 && n_to_find != count)
-    {
-      OCTAVE_QUIT;
-
-      if (nda(k) != T ())
-	{
-	  end_el = k;
-	  if (start_el == -1)
-	    start_el = k;
-	  count++;
-	}
-      k = k + direction;
-    }
-
-  // Reverse the range if we're looking backward.
-  if (direction == -1)
-    {
-      octave_idx_type tmp_el = start_el;
-      start_el = end_el;
-      end_el = tmp_el;
-    }
-  // Fix an off by one error.
-  end_el++;
-
-  // If the original argument was a row vector, force a row vector of
-  // the overall indices to be returned.  But see below for scalar
-  // case...
-
-  octave_idx_type result_nr = count;
-  octave_idx_type result_nc = 1;
-
-  bool column_vector_arg = false;
-  bool scalar_arg = false;
-
-  if (nda.ndims () == 2)
-    {
-      octave_idx_type nr = nda.rows ();
-      octave_idx_type nc = nda.columns ();
+  Array<octave_idx_type> idx;
+  if (n_to_find >= 0)
+    idx = nda.find (n_to_find, direction == -1);
+  else
+    idx = nda.find ();
 
-      if (nr == 1)
-	{
-	  result_nr = 1;
-	  result_nc = count;
-
-	  scalar_arg = (nc == 1);
-	}
-      else if (nc == 1)
-	column_vector_arg = true;
-    }
-
-  Matrix idx (result_nr, result_nc);
-
-  Matrix i_idx (result_nr, result_nc);
-  Matrix j_idx (result_nr, result_nc);
-
-  ArrayN<T> val (dim_vector (result_nr, result_nc));
-
-  if (count > 0)
-    {
-      count = 0;
-
-      octave_idx_type nr = nda.rows ();
-
-      octave_idx_type i = 0;
-
-      // Search for elements to return.  Only search the region where
-      // there are elements to be found using the count that we want
-      // to find.
+  // Fixup idx dimensions, for Matlab compatibility.
+  // find(zeros(0,0)) -> zeros(0,0)
+  // find(zeros(1,0)) -> zeros(1,0)
+  // find(zeros(0,1)) -> zeros(0,1)
+  // find(zeros(0,X)) -> zeros(0,1)
+  // find(zeros(1,1)) -> zeros(0,0) !!!! WHY?
+  // find(zeros(0,1,0)) -> zeros(0,0)
+  // find(zeros(0,1,0,1)) -> zeros(0,0) etc
+  // FIXME: I don't believe this is right. Matlab seems to violate its own docs
+  // here, because a scalar *is* a row vector.
 
-      // For compatibility, all N-d arrays are handled as if they are
-      // 2-d, with the number of columns equal to "prod (dims (2:end))".
-
-      for (k = start_el; k < end_el; k++)
-	{
-	  OCTAVE_QUIT;
-
-	  if (nda(k) != T ())
-	    {
-	      idx(count) = k + 1;
-
-	      octave_idx_type xr = k % nr;
-	      i_idx(count) = xr + 1;
-	      j_idx(count) = (k - xr) / nr + 1;
-
-	      val(count) = nda(k);
-
-	      count++;
-	    }
-
-	  i++;
-	}
-    }
-  else if (scalar_arg || (nda.rows () == 0 && ! column_vector_arg))
-    {
-      idx.resize (0, 0);
-
-      i_idx.resize (0, 0);
-      j_idx.resize (0, 0);
-
-      val.resize (dim_vector (0, 0));
-    }
+  if ((nda.numel () == 1 && idx.is_empty ())
+      || (nda.rows () == 0 && nda.dims ().numel (1) == 0))
+    idx = idx.reshape (dim_vector (0, 0));
+  else if (nda.rows () == 1 && nda.ndims () == 2)
+    idx = idx.reshape (dim_vector (1, idx.length ()));
 
   switch (nargout)
     {
     default:
     case 3:
-      retval(2) = val;
+      retval(2) = ArrayN<T> (nda.index (idx_vector (idx)));
       // Fall through!
 
     case 2:
-      retval(1) = j_idx;
-      retval(0) = i_idx;
-      break;
+      {
+        Array<octave_idx_type> jdx (idx.dims ());
+        octave_idx_type n = idx.length (), nr = nda.rows ();
+        for (octave_idx_type i = 0; i < n; i++)
+          {
+            jdx.xelem (i) = idx.xelem (i) / nr;
+            idx.xelem (i) %= nr;
+          }
+        retval(1) = NDArray (jdx, true);
+      }
+      // Fall through!
 
     case 1:
     case 0:
-      retval(0) = idx;
+      retval(0) = NDArray (idx, true);
       break;
     }
 
   return retval;
 }
 
-template octave_value_list find_nonzero_elem_idx (const Array<double>&, int,
-						  octave_idx_type, int);
-
-template octave_value_list find_nonzero_elem_idx (const Array<Complex>&, int,
-						  octave_idx_type, int);
+#define INSTANTIATE_FIND_ARRAY(T) \
+template octave_value_list find_nonzero_elem_idx (const Array<T>&, int, \
+						  octave_idx_type, int)
 
-template octave_value_list find_nonzero_elem_idx (const Array<float>&, int,
-						  octave_idx_type, int);
-
-template octave_value_list find_nonzero_elem_idx (const Array<FloatComplex>&,
-						  int, octave_idx_type, int);
+INSTANTIATE_FIND_ARRAY(double);
+INSTANTIATE_FIND_ARRAY(float);
+INSTANTIATE_FIND_ARRAY(Complex);
+INSTANTIATE_FIND_ARRAY(FloatComplex);
+INSTANTIATE_FIND_ARRAY(bool);
+INSTANTIATE_FIND_ARRAY(octave_int8);
+INSTANTIATE_FIND_ARRAY(octave_int16);
+INSTANTIATE_FIND_ARRAY(octave_int32);
+INSTANTIATE_FIND_ARRAY(octave_int64);
+INSTANTIATE_FIND_ARRAY(octave_uint8);
+INSTANTIATE_FIND_ARRAY(octave_uint16);
+INSTANTIATE_FIND_ARRAY(octave_uint32);
+INSTANTIATE_FIND_ARRAY(octave_uint64);
 
 template <typename T>
 octave_value_list
@@ -342,6 +262,9 @@
 template octave_value_list find_nonzero_elem_idx (const Sparse<Complex>&, int,
 						  octave_idx_type, int);
 
+template octave_value_list find_nonzero_elem_idx (const Sparse<bool>&, int,
+						  octave_idx_type, int);
+
 octave_value_list
 find_nonzero_elem_idx (const PermMatrix& v, int nargout, 
 		       octave_idx_type n_to_find, int direction)
@@ -561,7 +484,51 @@
 
   octave_value arg = args(0);
 
-  if (arg.is_sparse_type ())
+  if (arg.is_bool_type ())
+    {
+      if (arg.is_sparse_type ())
+        {
+	  SparseBoolMatrix v = arg.sparse_bool_matrix_value ();
+
+	  if (! error_state)
+	    retval = find_nonzero_elem_idx (v, nargout, 
+					    n_to_find, direction);
+        }
+      else
+        {
+          boolNDArray v = arg.bool_array_value ();
+
+	  if (! error_state)
+	    retval = find_nonzero_elem_idx (v, nargout, 
+					    n_to_find, direction);
+        }
+    }
+  else if (arg.is_integer_type ())
+    {
+#define DO_INT_BRANCH(INTT) \
+      else if (arg.is_ ## INTT ## _type ()) \
+        { \
+          INTT ## NDArray v = arg.INTT ## _array_value (); \
+          \
+	  if (! error_state) \
+	    retval = find_nonzero_elem_idx (v, nargout, \
+					    n_to_find, direction);\
+        }
+
+      if (false)
+        ;
+      DO_INT_BRANCH (int8)
+      DO_INT_BRANCH (int16)
+      DO_INT_BRANCH (int32)
+      DO_INT_BRANCH (int64)
+      DO_INT_BRANCH (uint8)
+      DO_INT_BRANCH (uint16)
+      DO_INT_BRANCH (uint32)
+      DO_INT_BRANCH (uint64)
+      else
+        panic_impossible ();
+    }
+  else if (arg.is_sparse_type ())
     {
       if (arg.is_real_type ())
 	{
--- a/src/version.h	Wed Mar 25 18:13:08 2009 -0700
+++ b/src/version.h	Thu Mar 26 09:47:35 2009 -0700
@@ -25,11 +25,11 @@
 #if !defined (octave_version_h)
 #define octave_version_h 1
 
-#define OCTAVE_VERSION "3.1.54+"
+#define OCTAVE_VERSION "3.1.55"
 
-#define OCTAVE_API_VERSION "api-v36+"
+#define OCTAVE_API_VERSION "api-v37"
 
-#define OCTAVE_RELEASE_DATE "2009-03-07"
+#define OCTAVE_RELEASE_DATE "2009-03-25"
 
 #define OCTAVE_COPYRIGHT "Copyright (C) 2009 John W. Eaton and others."