Mercurial > hg > octave-lojdl
changeset 9647:54f45f883a53
optimize & extend randperm
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Wed, 16 Sep 2009 13:41:49 +0200 |
parents | b6f5a59a66d7 |
children | 11844593875a |
files | liboctave/ChangeLog liboctave/oct-rand.cc scripts/ChangeLog scripts/general/Makefile.in scripts/general/randperm.m src/ChangeLog src/DLD-FUNCTIONS/rand.cc |
diffstat | 7 files changed, 103 insertions(+), 44 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,8 @@ +2009-09-16 Jaroslav Hajek <highegg@gmail.com> + + * oct-rand.cc (octave_rand::do_matrix, do_nd_array, do_vector): + Use Array::clear rather than Array::resize. + 2009-09-06 Jaroslav Hajek <highegg@gmail.com> * dColVector.h (operator *(const Matrix&, const ColumnVector)):
--- a/liboctave/oct-rand.cc +++ b/liboctave/oct-rand.cc @@ -390,7 +390,7 @@ if (n >= 0 && m >= 0) { - retval.resize (n, m); + retval.clear (n, m); if (n > 0 && m > 0) fill (retval.capacity(), retval.fortran_vec(), a); @@ -408,7 +408,7 @@ if (! dims.all_zero ()) { - retval.resize (dims); + retval.clear (dims); fill (retval.capacity(), retval.fortran_vec(), a); } @@ -423,7 +423,7 @@ if (n > 0) { - retval.resize (n); + retval.clear (n); fill (retval.capacity (), retval.fortran_vec (), a); }
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2009-09-16 Jaroslav Hajek <highegg@gmail.com> + + * general/randperm.m: Remove. + * general/Makefile.in: Update. + 2009-09-15 John W. Eaton <jwe@octave.org> * confiugre.ac: Rename from configure.in
--- a/scripts/general/Makefile.in +++ b/scripts/general/Makefile.in @@ -43,7 +43,7 @@ isequalwithequalnans.m isscalar.m issquare.m issymmetric.m \ isvector.m loadobj.m logical.m logspace.m mod.m nargchk.m \ nargoutchk.m nextpow2.m nthroot.m num2str.m perror.m pol2cart.m \ - polyarea.m postpad.m prepad.m quadgk.m quadl.m quadv.m randperm.m rat.m \ + polyarea.m postpad.m prepad.m quadgk.m quadl.m quadv.m rat.m \ rem.m repmat.m rot90.m rotdim.m runlength.m saveobj.m shift.m shiftdim.m \ sortrows.m sph2cart.m strerror.m structfun.m subsindex.m \ triplequad.m trapz.m tril.m triu.m
deleted file mode 100644 --- a/scripts/general/randperm.m +++ /dev/null @@ -1,40 +0,0 @@ -## Copyright (C) 1998, 1999, 2000, 2002, 2005, 2006, 2007 John W. Eaton -## -## This file is part of Octave. -## -## 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 3 of the License, 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 Octave; see the file COPYING. If not, see -## <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn {Function File} {} randperm (@var{n}) -## Return a row vector containing a random permutation of the -## integers from 1 to @var{n}. -## @end deftypefn - -## Author: "James R. Van Zandt" <jrv@vanzandt.mv.com> -## Adapted-By: jwe - -function retval = randperm (n) - - if (nargin == 1 && isscalar (n) && floor (n) == n) - if (n >= 0) - [junk, retval] = sort (rand (1, n)); - else - error ("randperm: argument must be non-negative"); - endif - else - print_usage (); - endif - -endfunction
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,7 @@ +2009-09-16 Jaroslav Hajek <highegg@gmail.com> + + * DLD-FUNCTIONS/rand.cc (Frandperm): New function. + 2009-09-15 Jaroslav Hajek <highegg@gmail.com> * pt-misc.cc (tree_parameter_list::convert_to_const_vector): Pass
--- a/src/DLD-FUNCTIONS/rand.cc +++ b/src/DLD-FUNCTIONS/rand.cc @@ -2,6 +2,7 @@ Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2003, 2005, 2006, 2007, 2008, 2009 John W. Eaton +Copyright (C) 2009 VZLU Prague This file is part of Octave. @@ -40,6 +41,7 @@ #include "oct-obj.h" #include "unwind-prot.h" #include "utils.h" +#include "ov-re-mat.h" /* %!shared __random_statistical_tests__ @@ -1017,6 +1019,89 @@ %! endif */ +DEFUN_DLD (randperm, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {} randperm (@var{n})\n\ +@deftypefnx {Loadable Function} {} randperm (@var{n}, @var{m})\n\ +Return a row vector containing a random permutation of @code{1:@var{n}}.\n\ +If @var{m} is supplied, return @var{m} permutations,\n\ +one in each row of a NxM matrix. The complexity is O(M*N) in both time and\n\ +memory. The randomization is performed using rand().\n\ +All permutations are equally likely.\n\ +@seealso{perms}\n\ +@end deftypefn") +{ + int nargin = args.length (); + octave_value retval; + + if (nargin == 1 || nargin == 2) + { + octave_idx_type n, m; + + if (nargin == 2) + m = args(1).idx_type_value (true); + else + m = 1; + + n = args(0).idx_type_value (true); + + if (m < 0 || n < 0) + error ("randperm: m and n must be non-negative"); + + if (! error_state) + { + // Generate random numbers. + NDArray r = octave_rand::nd_array (dim_vector (m, n)); + + // Create transposed to allow faster access. + Array<octave_idx_type> idx (dim_vector (n, m)); + + double *rvec = r.fortran_vec (); + + octave_idx_type *ivec = idx.fortran_vec (); + + // Perform the Knuth shuffle. + for (octave_idx_type j = 0; j < m; j++) + { + for (octave_idx_type i = 0; i < n; i++) + ivec[i] = i; + + for (octave_idx_type i = 0; i < n; i++) + { + octave_idx_type k = i + floor (rvec[i] * (n - i)); + std::swap (ivec[i], ivec[k]); + } + + ivec += n; + rvec += n; + } + + // Transpose. + idx = idx.transpose (); + + // Re-fetch the pointers. + ivec = idx.fortran_vec (); + rvec = r.fortran_vec (); + + // Convert to doubles, reusing r. + for (octave_idx_type i = 0, l = m*n; i < l; i++) + rvec[i] = ivec[i] + 1; + + // Now create an array object with a cached idx_vector. + retval = new octave_matrix (r, idx_vector (idx)); + } + } + else + print_usage (); + + return retval; +} + +/* +%!assert(sort(randperm(20)),1:20) +%!assert(sort(randperm(20,50),2),repmat(1:20,50,1)) +*/ + /* ;;; Local Variables: *** ;;; mode: C++ ***