# HG changeset patch # User Jaroslav Hajek # Date 1253101309 -7200 # Node ID 54f45f883a53a7805039d3906ad461856d965db5 # Parent b6f5a59a66d73144323a3578e266cbf23eeab5a8 optimize & extend randperm diff -r b6f5a59a66d7 -r 54f45f883a53 liboctave/ChangeLog --- a/liboctave/ChangeLog Tue Sep 15 21:19:15 2009 +0200 +++ b/liboctave/ChangeLog Wed Sep 16 13:41:49 2009 +0200 @@ -1,3 +1,8 @@ +2009-09-16 Jaroslav Hajek + + * oct-rand.cc (octave_rand::do_matrix, do_nd_array, do_vector): + Use Array::clear rather than Array::resize. + 2009-09-06 Jaroslav Hajek * dColVector.h (operator *(const Matrix&, const ColumnVector)): diff -r b6f5a59a66d7 -r 54f45f883a53 liboctave/oct-rand.cc --- a/liboctave/oct-rand.cc Tue Sep 15 21:19:15 2009 +0200 +++ b/liboctave/oct-rand.cc Wed Sep 16 13:41:49 2009 +0200 @@ -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); } diff -r b6f5a59a66d7 -r 54f45f883a53 scripts/ChangeLog --- a/scripts/ChangeLog Tue Sep 15 21:19:15 2009 +0200 +++ b/scripts/ChangeLog Wed Sep 16 13:41:49 2009 +0200 @@ -1,3 +1,8 @@ +2009-09-16 Jaroslav Hajek + + * general/randperm.m: Remove. + * general/Makefile.in: Update. + 2009-09-15 John W. Eaton * confiugre.ac: Rename from configure.in diff -r b6f5a59a66d7 -r 54f45f883a53 scripts/general/Makefile.in --- a/scripts/general/Makefile.in Tue Sep 15 21:19:15 2009 +0200 +++ b/scripts/general/Makefile.in Wed Sep 16 13:41:49 2009 +0200 @@ -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 diff -r b6f5a59a66d7 -r 54f45f883a53 scripts/general/randperm.m --- a/scripts/general/randperm.m Tue Sep 15 21:19:15 2009 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -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 -## . - -## -*- 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" -## 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 diff -r b6f5a59a66d7 -r 54f45f883a53 src/ChangeLog --- a/src/ChangeLog Tue Sep 15 21:19:15 2009 +0200 +++ b/src/ChangeLog Wed Sep 16 13:41:49 2009 +0200 @@ -1,3 +1,7 @@ +2009-09-16 Jaroslav Hajek + + * DLD-FUNCTIONS/rand.cc (Frandperm): New function. + 2009-09-15 Jaroslav Hajek * pt-misc.cc (tree_parameter_list::convert_to_const_vector): Pass diff -r b6f5a59a66d7 -r 54f45f883a53 src/DLD-FUNCTIONS/rand.cc --- a/src/DLD-FUNCTIONS/rand.cc Tue Sep 15 21:19:15 2009 +0200 +++ b/src/DLD-FUNCTIONS/rand.cc Wed Sep 16 13:41:49 2009 +0200 @@ -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 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++ ***