Mercurial > octave-nkf

changeset 9647:54f45f883a53

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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	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  <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	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);
     }
--- 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  <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	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
--- 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
-## <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	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  <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	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<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++ ***