# HG changeset patch
# User jwe
# Date 1068571542 0
# Node ID 77566be8b9e95c0c55045939bb9bcf0c6538ba50
# Parent  a97b498e1b321d37613e3e48b54a7d275693ed1b
[project @ 2003-11-11 17:25:42 by jwe]

diff --git a/liboctave/Array-util.cc b/liboctave/Array-util.cc
--- a/liboctave/Array-util.cc
+++ b/liboctave/Array-util.cc
@@ -510,6 +510,25 @@
   return retval;
 }
 
+Array<int>
+calc_permutated_idx (const Array<int>& old_idx, 
+		     const Array<int>& perm_vec, bool inv)
+{
+  int n_el = old_idx.length ();
+
+  Array<int> retval (n_el);
+
+  for (int i = 0; i < n_el; i++)
+    {
+      if (inv)
+	retval(perm_vec(i)-1) = old_idx(i);
+      else
+	retval(i) = old_idx(perm_vec(i)-1);
+    }
+
+  return retval;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
diff --git a/liboctave/Array-util.h b/liboctave/Array-util.h
--- a/liboctave/Array-util.h
+++ b/liboctave/Array-util.h
@@ -87,6 +87,8 @@
 				const dim_vector& dimensions,
 				int resize_ok);
 
+extern Array<int> calc_permutated_idx (const Array<int>& old_idx, 
+				       const Array<int>& perm_vec, bool inv);
 #endif
 
 /*
diff --git a/liboctave/Array.cc b/liboctave/Array.cc
--- a/liboctave/Array.cc
+++ b/liboctave/Array.cc
@@ -356,6 +356,71 @@
 }
 
 template <class T>
+Array<T>
+Array<T>::permute (const Array<int>& perm_vec, bool inv) const
+{
+  Array<T> retval;
+
+  dim_vector dv = dims ();
+  dim_vector dv_new;
+
+  int nd = dv.length ();
+
+  dv_new.resize (nd);
+
+  // Need this array to check for identical elements in permutation array.
+  Array<bool> checked (nd, false);
+
+  // Find dimension vector of permuted array.
+  for (int i = 0; i < nd; i++)
+    {
+      int perm_el = perm_vec.elem (i);
+
+      if (perm_el > dv.length () || perm_el < 1)
+	{
+	  (*current_liboctave_error_handler)
+	    ("permutation vector contains an invalid element");
+
+	  return retval;
+	}
+
+      if (checked.elem(perm_el - 1))
+	{
+	  (*current_liboctave_error_handler)
+	    ("PERM cannot contain identical elements");
+
+	  return retval;
+	}
+      else
+	checked.elem(perm_el - 1) = true;
+
+      dv_new (i) = dv (perm_el - 1);
+    }
+
+  retval.resize (dv_new);
+
+  // Index array to the original array.
+  Array<int> old_idx (nd, 0);
+
+  // Number of elements in Array (should be the same for
+  // both the permuted array and original array).
+  int n = retval.length ();
+
+  // Permute array.
+  for (int i = 0; i < n; i++)
+    {
+      // Get the idx of permuted array.
+      Array<int> new_idx = calc_permutated_idx (old_idx, perm_vec, inv);
+
+      retval.elem (new_idx) = elem (old_idx);
+
+      increment_index (old_idx, dv);
+    }
+
+  return retval;
+}
+
+template <class T>
 void
 Array<T>::resize_no_fill (int n)
 {
diff --git a/liboctave/Array.h b/liboctave/Array.h
--- a/liboctave/Array.h
+++ b/liboctave/Array.h
@@ -407,6 +407,10 @@
 
   Array<T> reshape (const dim_vector& new_dims) const;
 
+  Array<T> permute (const Array<int>& vec, bool inv = false) const;
+  Array<T> ipermute (const Array<int>& vec) const
+    { return permute (vec, true); }
+
   void resize_no_fill (int n);
   void resize_and_fill (int n, const T& val);
 
diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -2220,32 +2220,32 @@
   OCTAVE_QUIT;
 
   // construct balancing permutation vector
-  Array<int> ipermute (nc);
+  Array<int> iperm (nc);
   for (int i = 0; i < nc; i++)
-    ipermute(i) = i;  // initialize to identity permutation
+    iperm(i) = i;  // initialize to identity permutation
 
   // leading permutations in forward order
   for (int i = 0; i < (ilo-1); i++)
     {
       int swapidx = static_cast<int> (dpermute(i)) - 1;
-      int tmp = ipermute(i);
-      ipermute(i) = ipermute(swapidx);
-      ipermute(swapidx) = tmp;
+      int tmp = iperm(i);
+      iperm(i) = iperm(swapidx);
+      iperm(swapidx) = tmp;
     }
 
   // trailing permutations must be done in reverse order
   for (int i = nc - 1; i >= ihi; i--)
     {
       int swapidx = static_cast<int> (dpermute(i)) - 1;
-      int tmp = ipermute(i);
-      ipermute(i) = ipermute(swapidx);
-      ipermute(swapidx) = tmp;
+      int tmp = iperm(i);
+      iperm(i) = iperm(swapidx);
+      iperm(swapidx) = tmp;
     }
 
   // construct inverse balancing permutation vector
   Array<int> invpvec (nc);
   for (int i = 0; i < nc; i++)
-    invpvec(ipermute(i)) = i;     // Thanks to R. A. Lippert for this method
+    invpvec(iperm(i)) = i;     // Thanks to R. A. Lippert for this method
 
   OCTAVE_QUIT;
 
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,15 @@
+2003-11-11  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* Array.h (Array<T>::ipermute): New function.
+
+2003-11-11  Petter Risholm  <risholm@stud.ntnu.no>
+
+	* Array.cc (Array<T>::permute): New function.
+	* Array.h: Provide decl.
+
+	* Array-util.cc (calc_permutated_idx): New function.
+	* Array-util.h: Provide decl.
+
 2003-11-10  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* Array.cc (Array<T>::index2): Return value has orientation of
diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -1863,32 +1863,32 @@
   OCTAVE_QUIT;
 
   // construct balancing permutation vector
-  Array<int> ipermute (nc);
+  Array<int> iperm (nc);
   for (int i = 0; i < nc; i++)
-    ipermute(i) = i;  // identity permutation
+    iperm(i) = i;  // identity permutation
 
   // leading permutations in forward order
   for (int i = 0; i < (ilo-1); i++)
     {
       int swapidx = static_cast<int> (dpermute(i)) - 1;
-      int tmp = ipermute(i);
-      ipermute(i) = ipermute (swapidx);
-      ipermute(swapidx) = tmp;
+      int tmp = iperm(i);
+      iperm(i) = iperm (swapidx);
+      iperm(swapidx) = tmp;
     }
 
   // trailing permutations must be done in reverse order
   for (int i = nc - 1; i >= ihi; i--)
     {
       int swapidx = static_cast<int> (dpermute(i)) - 1;
-      int tmp = ipermute(i);
-      ipermute(i) = ipermute(swapidx);
-      ipermute(swapidx) = tmp;
+      int tmp = iperm(i);
+      iperm(i) = iperm(swapidx);
+      iperm(swapidx) = tmp;
     }
 
   // construct inverse balancing permutation vector
   Array<int> invpvec (nc);
   for (int i = 0; i < nc; i++)
-    invpvec(ipermute(i)) = i;     // Thanks to R. A. Lippert for this method
+    invpvec(iperm(i)) = i;     // Thanks to R. A. Lippert for this method
 
   OCTAVE_QUIT;
  
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,19 @@
+2003-11-11  John W. Eaton  <jwe@bevo.che.wisc.edu>
+
+	* data.cc (do_permute): New function.
+	(Fpermute, Fipermute): Use it.
+
+	* ov-base.cc (octave_base_value::permute): New function.
+	* ov-base.h: Provide decl.
+
+	* ov.h (octave_value::ipermute): New function.
+
+2003-11-11  Petter Risholm  <risholm@stud.ntnu.no>
+
+	* data.cc (Fpermute, Fipermute): New functions.
+	* ov.h (octave_value::permute): New function.
+	* ov-base-mat.h (octave_base_matrix::permute): New function.
+
 2003-11-10  John W. Eaton  <jwe@bevo.che.wisc.edu>
 
 	* oct-obj.h (octave_value_list): Internal representation is now
diff --git a/src/data.cc b/src/data.cc
--- a/src/data.cc
+++ b/src/data.cc
@@ -665,6 +665,55 @@
   DATA_REDUCTION (prod);
 }
 
+static octave_value
+do_permute (const octave_value_list& args, bool inv, const std::string& fname)
+{
+  octave_value retval;
+
+  if (args.length () == 2 && args(1).length () == args(0).dims ().length ())
+    {
+      Array<int> vec = args(1).int_vector_value ();
+
+      octave_value ret = args(0).permute (vec, inv);
+
+      if (! error_state)
+	retval = ret;
+    }
+  else
+    print_usage (fname);
+
+  return retval;
+}
+
+DEFUN (permute, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Built-in Function} {} permute (@var{a}, @var{perm})\n\
+Return the generalized transpose for an N-d array object @var{a}.\n\
+The permutation vector @var{perm} must contain the elements\n\
+@code{1:ndims(a)} (in any order, but each element must appear just once).\n\
+\n\
+@end deftypefn\n\
+@seealso{ipermute}")
+{
+  return do_permute (args, false, "permute");
+}
+
+DEFUN (ipermute, args, ,
+  "-*- texinfo -*-\n\
+@deftypefn {Built-in Function} {} ipermute (@var{a}, @var{iperm})\n\
+The inverse of the @code{permute} function.  The expression\n\
+\n\
+@example\n\
+ipermute (permute (a, perm), perm)\n\
+@end example\n\
+returns the original array @var{a}.\n\
+\n\
+@end deftypefn\n\
+@seealso{permute}")
+{
+  return do_permute (args, true, "ipermute");
+}
+
 DEFUN (length, args, ,
   "-*- texinfo -*-\n\
 @deftypefn {Built-in Function} {} length (@var{a})\n\
diff --git a/src/ov-base-mat.h b/src/ov-base-mat.h
--- a/src/ov-base-mat.h
+++ b/src/ov-base-mat.h
@@ -95,6 +95,9 @@
   octave_value reshape (const dim_vector& new_dims) const
     { return MT (matrix.reshape (new_dims)); }
 
+  octave_value permute (const Array<int>& vec, bool inv = false) const
+    { return MT (matrix.permute (vec, inv)); }
+
   octave_value all (int dim = 0) const { return matrix.all (dim); }
   octave_value any (int dim = 0) const { return matrix.any (dim); }
 
diff --git a/src/ov-base.cc b/src/ov-base.cc
--- a/src/ov-base.cc
+++ b/src/ov-base.cc
@@ -181,6 +181,13 @@
 }
 
 octave_value
+octave_base_value::permute (const Array<int>&, bool) const
+{
+  gripe_wrong_type_arg ("octave_base_value::permute ()", type_name ());
+  return octave_value ();
+}
+
+octave_value
 octave_base_value::convert_to_str_internal (bool, bool) const
 {
   gripe_wrong_type_arg ("octave_base_value::convert_to_str_internal ()",
diff --git a/src/ov-base.h b/src/ov-base.h
--- a/src/ov-base.h
+++ b/src/ov-base.h
@@ -98,6 +98,8 @@
 
   octave_value reshape (const dim_vector&) const;
 
+  octave_value permute (const Array<int>& vec, bool = false) const;
+
   bool is_defined (void) const { return false; }
 
   bool is_cell (void) const { return false; }
diff --git a/src/ov.h b/src/ov.h
--- a/src/ov.h
+++ b/src/ov.h
@@ -332,6 +332,12 @@
   virtual octave_value reshape (const dim_vector& dv) const
     { return rep->reshape (dv); }
 
+  virtual octave_value permute (const Array<int>& vec, bool inv = false) const
+    { return rep->permute (vec, inv); }
+
+  octave_value ipermute (const Array<int>& vec) const
+    { return rep->permute (vec, true); }
+
   // Does this constant have a type?  Both of these are provided since
   // it is sometimes more natural to write is_undefined() instead of
   // ! is_defined().