octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: find first and last


From: David Bateman
Subject: Re: find first and last
Date: Wed, 27 Sep 2006 16:21:30 +0200
User-agent: Thunderbird 1.5.0.7 (X11/20060921)

John W. Eaton wrote:
> On 27-Sep-2006, David Bateman wrote:
>
> | John W. Eaton wrote:
> | > Oh, it would also be nice if spfind also handled the N and
> | > "last"/"first" arguments.  So let me register that as a feature
> | > request.
> | > 
> | > jwe
> | > 
> | > 
> | 
> | What about something like the attached.
>
> If this provides compatibility with Matlab's find function for sparse
> matrices, then I'd say check it in as a replacement for the current
> spfind.  I see that the current version may also return NR and NC.  Is
> there any need to preserve that behavior?
>
> jwe
>
>   
Ok, please consider the attached patch.

D.

-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

2006-09-27  David Bateman  <address@hidden>

        * DLD-FUNCTIONS/sparse.cc (spfind, sparse_find): Delete.
        * DLD-FUNCTIONS/spfind.cc: New file implementating compatible
        sparse find function.
        * Makefile.in (DLD_XSRC): Add spfind.cc.
        
*** ./src/DLD-FUNCTIONS/spfind.cc.orig3 2006-09-27 16:14:39.662035119 +0200
--- ./src/DLD-FUNCTIONS/spfind.cc       2006-09-27 16:20:07.899692040 +0200
***************
*** 0 ****
--- 1,288 ----
+ /*
+ 
+ Copyright (C) 2006 David Bateman
+ 
+ Octave is free software; you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by the
+ Free Software Foundation; either version 2, or (at your option) any
+ later version.
+ 
+ Octave is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+ for more details.
+ 
+ You should have received a copy of the GNU General Public License
+ along with this program; see the file COPYING.  If not, write to the
+ Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ Boston, MA 02110-1301, USA.
+ 
+ */
+ 
+ #ifdef HAVE_CONFIG_H
+ #include <config.h>
+ #endif
+ 
+ #include <cstdlib>
+ #include <string>
+ 
+ #include "variables.h"
+ #include "utils.h"
+ #include "pager.h"
+ #include "defun-dld.h"
+ #include "gripes.h"
+ #include "quit.h"
+ 
+ #include "ov-re-sparse.h"
+ #include "ov-cx-sparse.h"
+ 
+ template <typename T, typename M>
+ octave_value_list
+ sparse_find_non_zero_elem_idx (const T& v, M& val, int nargout, 
+                             octave_idx_type n_to_find, int direction)
+ {
+   octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ());
+ 
+ 
+   octave_idx_type nc = v.cols();
+   octave_idx_type nr = v.rows();
+   octave_idx_type nz = v.nnz();
+ 
+   // Search in the default range.
+   octave_idx_type start_nc = -1;
+   octave_idx_type end_nc = -1;
+   octave_idx_type count;
+  
+   // Search for the range to search
+   if (n_to_find < 0)
+     {
+       start_nc = 0;
+       end_nc = nc;
+       n_to_find = nz;
+       count = nz;
+     }
+   else if (direction > 0)
+     {
+       for (octave_idx_type j = 0; j < nc; j++)
+       {
+         OCTAVE_QUIT;
+         if (v.cidx(j) == 0 && v.cidx(j+1) != 0)
+           start_nc = j;
+         if (v.cidx(j+1) >= n_to_find)
+           {
+             end_nc = j + 1;
+             break;
+           }
+       }
+     }
+   else
+     {
+       for (octave_idx_type j = nc; j > 0; j--)
+       {
+         OCTAVE_QUIT;
+         if (v.cidx(j) == nz && v.cidx(j-1) != nz)
+           end_nc = j;
+         if (nz - v.cidx(j-1) >= n_to_find)
+           {
+             start_nc = j - 1;
+             break;
+           }
+       }
+     }
+ 
+   count = (n_to_find > v.cidx(end_nc) - v.cidx(start_nc) ? 
+          v.cidx(end_nc) - v.cidx(start_nc) : n_to_find);
+ 
+   // 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 scalar_arg = false;
+ 
+   if (v.rows () == 1)
+     {
+       result_nr = 1;
+       result_nc = count;
+ 
+       scalar_arg = (v.columns () == 1);
+     }
+ 
+   Matrix idx (result_nr, result_nc);
+ 
+   Matrix i_idx (result_nr, result_nc);
+   Matrix j_idx (result_nr, result_nc);
+ 
+   val = M(result_nr, result_nc);
+ 
+   if (count > 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.
+       for (octave_idx_type j = start_nc, cx = 0; j < end_nc; j++) 
+       for (octave_idx_type i = v.cidx(j); i < v.cidx(j+1); i++ ) 
+         {
+           OCTAVE_QUIT;
+           if (direction < 0 && i < nz - count)
+             continue;
+           i_idx (cx) = static_cast<double> (v.ridx(i) + 1);
+           j_idx (cx) = static_cast<double> (j + 1);
+           idx (cx) = j * nc + v.ridx(i) + 1; 
+           val (cx) = v.data(i);
+           cx++;
+           if (cx == count)
+             break;
+         }
+     }
+   else if (scalar_arg)
+     {
+       idx.resize (0, 0);
+ 
+       i_idx.resize (0, 0);
+       j_idx.resize (0, 0);
+ 
+       val.resize (0, 0);
+     }
+ 
+   switch (nargout)
+     {
+     case 0:
+     case 1:
+       retval(0) = idx;
+       break;
+ 
+     case 5:
+       retval(4) = nc;
+       // Fall through
+ 
+     case 4:
+       retval(3) = nr;
+       // Fall through
+ 
+     case 3:
+       retval(2) = val;
+       // Fall through!
+ 
+     case 2:
+       retval(1) = j_idx;
+       retval(0) = i_idx;
+       break;
+ 
+     default:
+       panic_impossible ();
+       break;
+     }
+ 
+   return retval;
+ }
+ 
+ template octave_value_list sparse_find_non_zero_elem_idx 
+   (const SparseMatrix&, Matrix&, int, octave_idx_type, int);
+ 
+ template octave_value_list sparse_find_non_zero_elem_idx 
+   (const SparseComplexMatrix&, ComplexMatrix&, int, octave_idx_type, int);
+ 
+ 
+ // PKG_ADD: dispatch ("find", "spfind", "sparse matrix");
+ // PKG_ADD: dispatch ("find", "spfind", "sparse complex matrix");
+ // PKG_ADD: dispatch ("find", "spfind", "sparse bool matrix");
+ DEFUN_DLD (spfind, args, nargout,
+     "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} {} spfind (@var{x})\n\
+ @deftypefnx {Loadable Function} {} spfind (@var{x}, @var{n})\n\
+ @deftypefnx {Loadable Function} {} spfind (@var{x}, @var{n}, 
@var{direction})\n\
+ @deftypefnx {Loadable Function} address@hidden, @var{j}, @var{v}} spfind 
(@dots{})\n\
+ \n\
+ A sparse version of the @code{find} function. Please see the @code{find}\n\
+ for details of its use.\n\
+ \n\
+ Note that this function is particularly useful for sparse matrices, as\n\
+ it extracts the non-zero elements as vectors, which can then be used to\n\
+ create the original matrix. For example,\n\
+ \n\
+ @example\n\
+ @group\n\
+ sz = size(a);\n\
+ [i, j, v] = spfind (a);\n\
+ b = sparse(i, j, v, sz(1), sz(2));\n\
+ @end group\n\
+ @end example\n\
+ @seealso{sparse}\n\
+ @end deftypefn")
+ {
+    octave_value_list retval;
+    int nargin = args.length ();
+ 
+   if (nargin > 3 || nargin < 1)
+     {
+       print_usage ();
+       return retval;
+     }
+ 
+   // Setup the default options.
+   octave_idx_type n_to_find = -1;
+   if (nargin > 1)
+     {
+       n_to_find = args(1).int_value ();
+       if (error_state)
+       {
+         error ("find: expecting second argument to be an integer");
+         return retval;
+       }
+     }
+ 
+   // Direction to do the searching (1 == forward, -1 == reverse).
+   int direction = 1;
+   if (nargin > 2)
+     {
+       direction = 0;
+ 
+       std::string s_arg = args(2).string_value ();
+ 
+       if (! error_state)
+       {
+         if (s_arg == "first")
+           direction = 1;
+         else if (s_arg == "last")
+           direction = -1;
+       }
+ 
+       if (direction == 0)
+       {
+         error ("find: expecting third argument to be \"first\" or \"last\"");
+         return retval;
+       }
+     }
+ 
+    octave_value arg = args(0);
+ 
+    if (arg.is_real_type ())
+      {
+        SparseMatrix v = arg.sparse_matrix_value ();
+        Matrix val;
+ 
+        if (! error_state)
+        retval = sparse_find_non_zero_elem_idx  (v, val, nargout, n_to_find, 
direction);
+      }
+    else if (arg.is_complex_type ())
+      {
+        SparseComplexMatrix v = arg.sparse_complex_matrix_value ();
+        ComplexMatrix val;
+ 
+        if (! error_state)
+        retval = sparse_find_non_zero_elem_idx  (v, val, nargout, n_to_find, 
direction);
+      }
+    else 
+      gripe_wrong_type_arg ("spfind", arg);
+ 
+    return retval;
+ }
+ 
+ /*
+ ;;; Local Variables: ***
+ ;;; mode: C++ ***
+ ;;; End: ***
+ */
*** ./src/DLD-FUNCTIONS/sparse.cc.orig3 2006-09-27 16:13:17.461417926 +0200
--- ./src/DLD-FUNCTIONS/sparse.cc       2006-09-27 16:12:57.227158396 +0200
***************
*** 399,591 ****
    return retval;
  }
  
- static octave_value_list
- sparse_find (const SparseMatrix& v)
- {
-   octave_value_list retval;
-   octave_idx_type nnz = v.nnz ();
-   dim_vector dv = v.dims ();
-   octave_idx_type nr = dv(0);
-   octave_idx_type nc = dv (1);
- 
-   ColumnVector I (nnz), J (nnz);
-   ColumnVector S (nnz);
- 
-   for (octave_idx_type i = 0, cx = 0; i < nc; i++) 
-     {
-       OCTAVE_QUIT;
-       for (octave_idx_type j = v.cidx(i); j < v.cidx(i+1); j++ ) 
-       {
-         I (cx) = static_cast<double> (v.ridx(j) + 1);
-         J (cx) = static_cast<double> (i + 1);
-         S (cx) = v.data(j);
-         cx++;
-       }
-     }
- 
-   if (dv(0) == 1)
-     {
-       retval(0) = I.transpose ();
-       retval(1) = J.transpose ();
-       retval(2) = S.transpose ();
-     }
-   else
-     {
-       retval(0) = I;
-       retval(1) = J;
-       retval(2) = S;
-     }
-   retval(3) = static_cast<double> (nr);
-   retval(4) = static_cast<double> (nc);
- 
-   return retval;
- }
- 
- static octave_value_list
- sparse_find (const SparseComplexMatrix& v)
- {
-   octave_value_list retval;
-   octave_idx_type nnz = v.nnz ();
-   dim_vector dv = v.dims ();
-   octave_idx_type nr = dv(0);
-   octave_idx_type nc = dv (1);
- 
-   ColumnVector I (nnz), J (nnz);
-   ComplexColumnVector S (nnz);
- 
-   for (octave_idx_type i = 0, cx = 0; i < nc; i++) 
-     {
-       OCTAVE_QUIT;
-       for (octave_idx_type j = v.cidx(i); j < v.cidx(i+1); j++ ) 
-       {
-         I (cx) = static_cast<double> (v.ridx(j) + 1);
-         J (cx) = static_cast<double> (i + 1);
-         S (cx) = v.data(j);
-         cx++;
-       }
-     }
- 
-   if (dv(0) == 1)
-     {
-       retval(0) = I.transpose ();
-       retval(1) = J.transpose ();
-       retval(2) = S.transpose ();
-     }
-   else
-     {
-       retval(0) = I;
-       retval(1) = J;
-       retval(2) = S;
-     }
-   retval(3) = static_cast<double> (nr);
-   retval(4) = static_cast<double> (nc);
- 
-   return retval;
- }
- 
- static octave_value_list
- sparse_find (const SparseBoolMatrix& v)
- {
-   octave_value_list retval;
-   octave_idx_type nnz = v.nnz ();
-   dim_vector dv = v.dims ();
-   octave_idx_type nr = dv(0);
-   octave_idx_type nc = dv (1);
- 
-   ColumnVector I (nnz), J (nnz);
-   ColumnVector S (nnz);
- 
-   for (octave_idx_type i = 0, cx = 0; i < nc; i++) 
-     {
-       OCTAVE_QUIT;
-       for (octave_idx_type j = v.cidx(i); j < v.cidx(i+1); j++ ) 
-       {
-         I (cx) = static_cast<double> (v.ridx(j) + 1);
-         J (cx) = static_cast<double> (i + 1);
-         S (cx) = static_cast<double> (v.data(j));
-         cx++;
-       }
-     }
- 
-   if (dv(0) == 1)
-     {
-       retval(0) = I.transpose ();
-       retval(1) = J.transpose ();
-       retval(2) = S.transpose ();
-     }
-   else
-     {
-       retval(0) = I;
-       retval(1) = J;
-       retval(2) = S;
-     }
-   retval(3) = static_cast<double> (nr);
-   retval(4) = static_cast<double> (nc);
- 
-   return retval;
- }
- 
- // PKG_ADD: dispatch ("find", "spfind", "sparse matrix");
- // PKG_ADD: dispatch ("find", "spfind", "sparse complex matrix");
- // PKG_ADD: dispatch ("find", "spfind", "sparse bool matrix");
- DEFUN_DLD (spfind, args, nargout ,
-     "-*- texinfo -*-\n\
- @deftypefn {Loadable Function} {[...] =} spfind (...)\n\
- SPFIND: a sparse version of the find operator\n\
- @enumerate\n\
-     @item\n\
- @var{x }= spfind( @var{a })\n\
-     @itemize @w\n\
- is analagous to @var{x}= find(@var{A}(:))@*\n\
- where @var{A}= full(@var{a})\n\
-     @end itemize\n\
-     @item\n\
- address@hidden,@var{j},@var{v},@var{nr},@var{nc}] = spfind( @var{a} )\n\
-     @itemize @w\n\
- returns column vectors @var{i},@var{j},@var{v} such address@hidden
- @var{a}= sparse(@var{i},@var{j},@var{v},@var{nr},@var{nc})\n\
-     @end itemize\n\
- @end enumerate\n\
- @seealso{sparse}\n\
- @end deftypefn")
- {
-    octave_value_list retval;
-    int nargin = args.length ();
- 
-    if (nargin != 1) 
-      {
-        print_usage ();
-        return retval;
-      }
-       
- 
-    octave_value arg = args(0);
- 
-    if (arg.is_sparse_type ())
-      {
-        if (arg.type_name () == "sparse matrix")
-        retval = sparse_find (args(0).sparse_matrix_value ());
-        else if (arg.type_name () == "sparse complex matrix" ) 
-        retval = sparse_find (args(0).sparse_complex_matrix_value ());
-        else if (arg.type_name () == "sparse bool matrix" ) 
-        retval = sparse_find (args(0).sparse_bool_matrix_value ());
-        else 
-        gripe_wrong_type_arg ("spfind", arg);
-      }
-    else
-      gripe_wrong_type_arg ("spfind", arg);
-      
-    if (nargout == 1 || nargout ==0 ) 
-      { 
-        // only find location as fortran index
-        octave_value_list tmp;
-        tmp(0) = retval(0) + (retval(1)-1)*retval(3);
-        retval = tmp;
-      }
- 
-    return retval;
- }
- 
  #define SPARSE_DIM_ARG_BODY(NAME, FUNC) \
      int nargin = args.length(); \
      octave_value retval; \
--- 399,404 ----
*** ./src/Makefile.in.orig3     2006-09-27 16:14:30.079227359 +0200
--- ./src/Makefile.in   2006-09-27 16:14:17.169527860 +0200
***************
*** 51,57 ****
        givens.cc hess.cc inv.cc kron.cc lpsolve.cc lsode.cc \
        lu.cc luinc.cc matrix_type.cc minmax.cc pinv.cc qr.cc \
        quad.cc qz.cc rand.cc regexp.cc schur.cc sort.cc sparse.cc \
!       spchol.cc spdet.cc spkron.cc splu.cc spparms.cc spqr.cc \
        sqrtm.cc svd.cc syl.cc time.cc \
        __gnuplot_raw__.l __glpk__.cc __pchip_deriv__.cc __qp__.cc
  
--- 51,57 ----
        givens.cc hess.cc inv.cc kron.cc lpsolve.cc lsode.cc \
        lu.cc luinc.cc matrix_type.cc minmax.cc pinv.cc qr.cc \
        quad.cc qz.cc rand.cc regexp.cc schur.cc sort.cc sparse.cc \
!       spchol.cc spdet.cc spfind.cc spkron.cc splu.cc spparms.cc spqr.cc \
        sqrtm.cc svd.cc syl.cc time.cc \
        __gnuplot_raw__.l __glpk__.cc __pchip_deriv__.cc __qp__.cc
  

reply via email to

[Prev in Thread] Current Thread [Next in Thread]