[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
- Re: find first and last, (continued)
- Re: find first and last, Bill Denney, 2006/09/26
- Re: find first and last, David Bateman, 2006/09/26
- Re: find first and last, Bill Denney, 2006/09/26
- Re: find first and last, John W. Eaton, 2006/09/26
- Re: find first and last, Bill Denney, 2006/09/26
- Re: find first and last, John W. Eaton, 2006/09/26
- Re: find first and last, John W. Eaton, 2006/09/26
- Re: find first and last, David Bateman, 2006/09/26
- Re: find first and last, John W. Eaton, 2006/09/27
- Re: find first and last, David Bateman, 2006/09/27
- Re: find first and last,
David Bateman <=
- Re: find first and last, John W. Eaton, 2006/09/27