Menu
▾
▴
matlisp-commit
|
From: Akshay S. <ak...@us...> - 2012-07-04 11:10:35
|
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "matlisp".
The branch, tensor has been updated
via d8e8b94a89920c6c741031b0a525fec2c62a9d2d (commit)
via 3727a088ffe014773472fd37a7d45346917a73a0 (commit)
from 695636685fd91ce1602b135d0c0e782ca06d47e7 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit d8e8b94a89920c6c741031b0a525fec2c62a9d2d
Author: Akshay Srinivasan <aks...@gm...>
Date: Wed Jul 4 16:36:01 2012 +0530
Cleaned up permutation.lisp
diff --git a/matlisp.asd b/matlisp.asd
index e300c13..cbdbae3 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -93,6 +93,8 @@
:depends-on ("standard-tensor"))
(:file "blas-helpers"
:depends-on ("standard-tensor"))
+ (:file "permutation"
+ :depends-on ("standard-tensor"))
;;
(:file "real-tensor"
:depends-on ("standard-tensor"))
diff --git a/src/conditions.lisp b/src/conditions.lisp
index ad16f1e..27986e0 100644
--- a/src/conditions.lisp
+++ b/src/conditions.lisp
@@ -13,78 +13,86 @@
(in-package :matlisp)
-;;---------------------------------------------------------------;;
-(define-condition generic-error (error)
- ((message :reader message :initarg :message :initform "")))
-
-(defmethod print-object ((c generic-error) stream)
- (format stream (message c)))
-
-;;---------------------------------------------------------------;;
-(define-condition invalid-type (generic-error)
+(defmacro defcondition (name (&rest parent-types) (&rest slot-specs) &body options)
+ "Like define-condition except that you can define
+ methods inside the condition definition with:
+ (:method {generic-function-name} {args*} &rest code*)
+"
+ (labels ((get-methods (opts mth rst)
+ (if (null opts) (values (reverse mth) (reverse rst))
+ (if (and (consp (car opts)) (eq (caar opts) :method))
+ (get-methods (cdr opts) (cons (car opts) mth) rst)
+ (get-methods (cdr opts) mth (cons (car opts) rst))))))
+ (multiple-value-bind (methods rest) (get-methods options nil nil)
+ `(progn
+ (define-condition ,name ,parent-types
+ ,slot-specs
+ ,@rest)
+ ,@(loop for mth in methods
+ collect `(defmethod ,@(cdr mth)))))))
+;;Generic conditions---------------------------------------------;;
+(defcondition generic-error (error)
+ ((message :reader message :initarg :message :initform ""))
+ (:method print-method ((c generic-error) stream)
+ (format stream (message c))))
+
+(defcondition invalid-type (generic-error)
((given-type :reader given :initarg :given)
(expected-type :reader expected :initarg :expected))
- (:documentation "Given an unexpected type."))
-
-(defmethod print-object ((c invalid-type) stream)
- (format stream "Given object of type ~A, expected ~A.~%" (given c) (expected c))
- (call-next-method))
+ (:documentation "Given an unexpected type.")
+ (:method print-object ((c invalid-type) stream)
+ (format stream "Given object of type ~A, expected ~A.~%" (given c) (expected c))
+ (call-next-method)))
-;;---------------------------------------------------------------;;
-(define-condition invalid-value (generic-error)
+(defcondition invalid-value (generic-error)
((given-value :reader given :initarg :given)
(expected-value :reader expected :initarg :expected))
- (:documentation "Given an unexpected value."))
+ (:documentation "Given an unexpected value.")
+ (:method print-object ((c invalid-value) stream)
+ (format stream "Given object ~A, expected ~A.~%" (given c) (expected c))
+ (call-next-method)))
-(defmethod print-object ((c invalid-value) stream)
- (format stream "Given object ~A, expected ~A.~%" (given c) (expected c))
- (call-next-method))
-
-;;---------------------------------------------------------------;;
-(define-condition unknown-token (generic-error)
+(defcondition unknown-token (generic-error)
((token :reader token :initarg :token))
- (:documentation "Given an unknown token."))
-
-(defmethod print-object ((c unknown-token) stream)
- (format stream "Given unknown token: ~A.~%" (token c))
- (call-next-method))
+ (:documentation "Given an unknown token.")
+ (:method print-object ((c unknown-token) stream)
+ (format stream "Given unknown token: ~A.~%" (token c))
+ (call-next-method)))
-;;---------------------------------------------------------------;;
-(define-condition coercion-error (generic-error)
+(defcondition coercion-error (generic-error)
((from :reader from :initarg :from)
(to :reader to :initarg :to))
- (:documentation "Cannot coerce one type into another."))
-
-(defmethod print-object ((c coercion-error) stream)
- (format stream "Cannot coerce ~a into ~a." (from c) (to c))
- (call-next-method))
+ (:documentation "Cannot coerce one type into another.")
+ (:method print-object ((c coercion-error) stream)
+ (format stream "Cannot coerce ~a into ~a." (from c) (to c))
+ (call-next-method)))
-;;---------------------------------------------------------------;;
-(define-condition matlisp-error (error)
+;;Tensor conditions----------------------------------------------;;
+(define-condition tensor-error (error)
;;Optional argument for error-handling.
((tensor :reader tensor :initarg :tensor)))
-(define-condition store-index-out-of-bounds (matlisp-error)
+(define-condition tensor-store-index-out-of-bounds (tensor-error)
((index :reader index :initarg :index)
(store-size :reader store-size :initarg :store-size))
(:documentation "An out of bounds index error for the one-dimensional store.")
(:report (lambda (c stream)
(format stream "Requested index ~A, but store is only of size ~A." (index c) (store-size c)))))
-(define-condition tensor-not-matrix (matlisp-error)
- ((tensor-rank :reader rank :initarg :rank))
- (:documentation "Given tensor is not a matrix.")
- (:report (lambda (c stream)
- (format stream "Given tensor with rank ~A, is not a matrix." (rank c)))))
-
-(define-condition insufficient-store (matlisp-error)
+(define-condition tensor-insufficient-store (tensor-error)
((store-size :reader store-size :initarg :store-size)
(max-idx :reader max-idx :initarg :max-idx))
(:documentation "Store is too small for the tensor with given dimensions.")
(:report (lambda (c stream)
(format stream "Store size is ~A, but maximum possible index is ~A." (store-size c) (max-idx c)))))
-(define-condition tensor-index-out-of-bounds (matlisp-error)
+(define-condition tensor-not-matrix (tensor-error)
+ ((tensor-rank :reader rank :initarg :rank))
+ (:documentation "Given tensor is not a matrix.")
+ (:report (lambda (c stream)
+ (format stream "Given tensor with rank ~A, is not a matrix." (rank c)))))
+
+(define-condition tensor-index-out-of-bounds (tensor-error)
((argument :reader argument :initarg :argument)
(index :reader index :initarg :index)
(argument-space-dimension :reader dimension :initarg :dimension))
@@ -92,48 +100,68 @@
(:report (lambda (c stream)
(format stream "~&Out of bounds for argument ~A: requested ~A, but dimension is only ~A." (argument c) (index c) (dimension c)))))
-(define-condition tensor-index-rank-mismatch (matlisp-error)
+(define-condition tensor-index-rank-mismatch (tensor-error)
((index-rank :reader index-rank :initarg :index-rank)
(rank :reader rank :initarg :rank))
(:documentation "Incorrect number of subscripts for the tensor.")
(:report (lambda (c stream)
(format stream "Index is of size ~A, whereas the tensor is of rank ~A." (index-rank c) (rank c)))))
-(define-condition tensor-invalid-head-value (matlisp-error)
+(define-condition tensor-invalid-head-value (tensor-error)
((head :reader head :initarg :head))
(:documentation "Incorrect value for the head of the tensor storage.")
(:report (lambda (c stream)
(format stream "Head of the store must be >= 0, initialized with ~A." (head c)))))
-(define-condition tensor-invalid-dimension-value (matlisp-error)
+(define-condition tensor-invalid-dimension-value (tensor-error)
((argument :reader argument :initarg :argument)
(argument-dimension :reader dimension :initarg :dimension))
(:documentation "Incorrect value for one of the dimensions of the tensor.")
(:report (lambda (c stream)
(format stream "Dimension of argument ~A must be > 0, initialized with ~A." (argument c) (dimension c)))))
-(define-condition tensor-invalid-stride-value (matlisp-error)
+(define-condition tensor-invalid-stride-value (tensor-error)
((argument :reader argument :initarg :argument)
(argument-stride :reader stride :initarg :stride))
(:documentation "Incorrect value for one of the strides of the tensor storage.")
(:report (lambda (c stream)
(format stream "Stride of argument ~A must be >= 0, initialized with ~A." (argument c) (stride c)))))
-(define-condition tensor-cannot-find-sub-class (matlisp-error)
+(define-condition tensor-cannot-find-sub-class (tensor-error)
((tensor-class :reader tensor-class :initarg :tensor-class))
(:documentation "Cannot find sub-class of the given tensor class")
(:report (lambda (c stream)
(format stream "Cannot find sub-class of the given tensor class: ~a." (tensor-class c)))))
-(define-condition tensor-cannot-find-optimization (matlisp-error)
+(define-condition tensor-cannot-find-optimization (tensor-error)
((tensor-class :reader tensor-class :initarg :tensor-class))
(:documentation "Cannot find optimization information for the given tensor class")
(:report (lambda (c stream)
(format stream "Cannot find optimization information for the given tensor class: ~a." (tensor-class c)))))
-(define-condition tensor-dimension-mismatch (matlisp-error)
+(define-condition tensor-dimension-mismatch (tensor-error)
()
(:documentation "The dimensions of the given tensors are not suitable for continuing with the operation.")
(:report (lambda (c stream)
(declare (ignore c))
(format stream "The dimensions of the given tensors are not suitable for continuing with the operation."))))
+
+;;Permutation conditions-----------------------------------------;;
+(define-condition permutation-error (error)
+ ((permutation :reader permutation :initarg :permutation)))
+
+(define-condition permutation-invalid-error (permutation-error)
+ ()
+ (:documentation "Object is not a permutation.")
+ (:report (lambda (c stream)
+ (declare (ignore c))
+ (format stream "Object is not a permutation."))))
+
+(define-condition permutation-permute-error (permutation-error)
+ ((sequence-length :reader seq-len :initarg :seq-len)
+ (group-rank :reader group-rank :initarg :group-rank))
+ (:documentation "Cannot permute sequence.")
+ (:report (lambda (c stream)
+ (format stream "Cannot permute sequence.
+sequence-length : ~a
+group-rank: ~a" (seq-len c) (group-rank c)))))
diff --git a/src/loopy-tests.lisp b/src/loopy-tests.lisp
index fbaafe0..8998588 100644
--- a/src/loopy-tests.lisp
+++ b/src/loopy-tests.lisp
@@ -44,7 +44,7 @@
do (setf (aref st-a of-a) (random 1d0)
(aref st-b of-b) (random 1d0)
(aref st-c of-c) 0d0)))
- (time
+ (time
(very-quickly
(mod-dotimes (idx (idxv n n n))
with (loop-order :row-major)
@@ -52,7 +52,7 @@
(of-a (idxv n 1 0))
(of-b (idxv 0 n 1))
(of-c (idxv n 0 1)))
- do (incf (aref st-c of-c) (* (aref st-a of-a) (aref st-b of-b)))))))))
+ do (incf (aref st-c of-c) (* (aref st-a of-a) (aref st-b of-b)))))))))
(defun test-mm-ddot (n)
(let* ((t-a (make-real-tensor-dims n n))
diff --git a/src/permutation.lisp b/src/permutation.lisp
new file mode 100644
index 0000000..40e0af2
--- /dev/null
+++ b/src/permutation.lisp
@@ -0,0 +1,317 @@
+(in-package :matlisp)
+
+(defun id-action-repr (n)
+ (declare (type fixnum n))
+ (let ((ret (allocate-index-store n)))
+ (declare (type (index-array *) ret))
+ (very-quickly
+ (loop
+ for i of-type index-type from 0 below n
+ do (setf (aref ret i) i)))
+ ret))
+
+(defun seqrnd (seq)
+ "Randomize the elements of a sequence. Destructive on SEQ."
+ (sort seq #'> :key #'(lambda (x) (random 1.0))))
+
+;;Class definitions----------------------------------------------;;
+(defclass permutation ()
+ ((representation :accessor repr
+ :initarg :repr)
+ (group-rank :accessor group-rank
+ :type index-type)))
+
+(defparameter +permutation-identity+
+ (let ((ret (make-instance 'permutation :repr :id)))
+ (setf (group-rank ret) 0)
+ ret))
+
+(defmethod print-object ((per permutation) stream)
+ (print-unreadable-object (per stream :type t)
+ (format stream "S_~a~%~a~%" (group-rank per) (repr per))))
+;;
+(defclass permutation-cycle (permutation)
+ ((representation :type cons)))
+
+(defun cycle-repr-p (perm)
+ "
+ Does a sorting operation to check for duplicate elements in
+ the cycle representation of a permutation.
+"
+ (if (not (typep perm '(index-array *))) nil
+ (locally
+ (declare (type (index-array *) perm))
+ (let ((len (length perm)))
+ (declare (type index-type len))
+ (if (<= len 1) nil
+ (let ((sort (very-quickly (sort (copy-seq perm) #'<))))
+ (declare (type (index-array *) sort))
+ (very-quickly
+ (loop for i of-type index-type from 1 below len
+ when (= (aref sort i) (aref sort (1- i)))
+ do (return nil)
+ finally (return t)))))))))
+
+(defmethod initialize-instance :after ((per permutation-cycle) &rest initargs)
+ (declare (ignore initargs))
+ (very-quickly
+ (loop
+ for cyc of-type (index-array *) in (repr per)
+ unless (cycle-repr-p cyc)
+ do (error 'permutation-invalid-error)
+ maximizing (idx-max cyc) into g-rnk of-type index-type
+ finally (setf (group-rank per) (the index-type (1+ g-rnk))))))
+
+(definline make-pcycle (&rest args)
+ (make-instance 'permutation-cycle :repr args))
+
+;;
+(defclass permutation-action (permutation)
+ ((representation :type (index-array *))))
+
+(defun action-repr-p (act)
+ "
+ Checks if ARR is a possible permutation vector. A permutation pi
+ is characterized by a vector containing the indices from 0,...,
+ @function{length}(@arg{perm})-1 in some order.
+"
+ (if (not (typep act '(index-array *))) nil
+ (locally
+ (declare (type (index-array *) act))
+ (let* ((len (length act))
+ (sort (very-quickly (sort (copy-seq act) #'<))))
+ (declare (type (index-array *) sort)
+ (type index-type len))
+ (very-quickly
+ (loop for i of-type index-type from 0 below len
+ unless (= (aref sort i) i)
+ do (return nil)
+ finally (return t)))))))
+
+(defmethod initialize-instance :after ((per permutation-action) &rest initargs)
+ (declare (ignore initargs))
+ (let ((act (repr per)))
+ (declare (type (index-array *) act))
+ (unless (action-repr-p act)
+ (error 'permutation-invalid-error))
+ (setf (group-rank per) (1+ (idx-max act)))))
+
+(definline make-paction (pact)
+ (make-instance 'permutation-action :repr pact))
+
+;;Conversions and validation-------------------------------------;;
+(defun action->cycle (act)
+ "
+ (action->cycle act)
+
+ This function obtains the canonical cycle representation
+ of a permutation. The first argument \"act\" is the action of the
+ permutation on the array #(0 1 2 3 ..): an object of the class
+ permutation-action.
+
+ \"Canonical\" may be a bit of an overstatement; this is the way
+ S_n was presented in Van der Waerden's book.
+"
+ (declare (type permutation-action act))
+ (mlet*
+ ((arr (repr act) :type (index-array *)))
+ (labels ((find-cycle (x0)
+ ;; This function obtains the cycle starting from x_0.
+ (declare (type index-type x0))
+ (if (= (aref arr x0) x0) (values 0 nil)
+ (very-quickly
+ (loop
+ for x of-type index-type = (aref arr x0) then (aref arr x)
+ and ret of-type cons = (list x0) then (cons x ret)
+ counting t into i of-type index-type
+ when (= x x0)
+ do (return (values i ret))))))
+ (cycle-walk (cyc ignore)
+ ;; Finds all cycles
+ (let ((x0 (find-if-not #'(lambda (x) (member x ignore)) arr)))
+ (if (null x0)
+ cyc
+ (multiple-value-bind (clen clst) (find-cycle x0)
+ (declare (type index-type clen)
+ (type list clst))
+ (cycle-walk
+ (if (= clen 0) cyc
+ (cons (make-array clen :element-type 'index-type :initial-contents clst) cyc))
+ (nconc ignore (if (= clen 0) (list x0) clst))))))))
+ (let ((cyc-lst (cycle-walk nil nil)))
+ (if (null cyc-lst)
+ +permutation-identity+
+ (make-instance 'permutation-cycle
+ :repr cyc-lst))))))
+
+(defun cycle->action (cyc)
+ "
+ (cycle->action cyc)
+
+ This function obtains the action representation of a permutation
+ from the cyclic one. The first argument \"cyc\" is the cyclic
+ representation of the permutation: an object of the class
+ permutation-cycle.
+"
+ (declare (type permutation-cycle cyc))
+ (let ((act-repr (id-action-repr (group-rank cyc)))
+ (cycs-repr (repr cyc)))
+ (declare (type (index-array *) act-repr))
+ (dolist (cyc cycs-repr)
+ (declare (type (index-array *) cyc))
+ (let ((xl (aref act-repr (aref cyc (1- (length cyc))))))
+ (very-quickly
+ (loop
+ for i of-type index-type downfrom (1- (length cyc)) to 0
+ do (setf (aref act-repr (aref cyc i))
+ (if (= i 0) xl
+ (aref act-repr (aref cyc (1- i)))))))))
+ (make-instance 'permutation-action :repr act-repr)))
+
+;;
+(defgeneric permute! (seq perm)
+ (:documentation "
+ (permute! seq perm)
+
+ Applies the permutation on the sequence.
+")
+ (:method ((seq sequence) (perm (eql +permutation-identity+)))
+ seq))
+
+(defmethod permute! ((seq sequence) (perm permutation-cycle))
+ (labels ((apply-cycle! (seq pcyc)
+ (declare (type (index-array *) pcyc))
+ (very-quickly
+ (let ((xl (aref seq (aref pcyc (1- (length pcyc))))))
+ (loop for i of-type index-type downfrom (1- (length pcyc)) to 0
+ do (setf (aref seq (aref pcyc i))
+ (if (= i 0) xl
+ (aref seq (aref pcyc (1- i))))))))))
+ (let ((len (length seq))
+ (glen (group-rank perm))
+ (cycs-lst (repr perm)))
+ (declare (type index-type len glen))
+ (if (< len glen) (error 'permutation-permute-error :seq-len len :group-rank glen)
+ (etypecase seq
+ (vector
+ (dolist (cyc cycs-lst seq)
+ (declare (type (index-array *) cyc))
+ (apply-cycle! seq cyc)))
+ (cons
+ (let ((cseq (make-array len :initial-contents seq)))
+ (declare (type (simple-vector *) cseq))
+ (dolist (cyc cycs-lst)
+ (declare (type (index-array *) cyc))
+ (apply-cycle! cseq cyc))
+ (mapl
+ (let ((i 0))
+ (declare (type fixnum i))
+ (lambda (x)
+ (when (< i glen)
+ (rplaca x (aref cseq i))
+ (incf i))))
+ seq))))))))
+
+(defmethod permute! ((seq sequence) (perm permutation-action))
+ (let ((len (length seq))
+ (glen (group-rank perm)))
+ (declare (type index-type len glen))
+ (if (< len glen) (error 'permutation-permute-error :seq-len len :group-rank glen)
+ (let ((cseq (make-array len :initial-contents seq))
+ (act (repr perm)))
+ (declare (type (simple-vector *) cseq)
+ (type (index-array *) act))
+ (etypecase seq
+ (vector
+ (very-quickly
+ (loop
+ for i of-type index-type from 0 below glen
+ do (unless (= i (aref act i))
+ (setf (aref seq i) (aref cseq (aref act i))))
+ finally (return seq))))
+ (cons
+ (mapl
+ (let ((i 0))
+ (declare (type fixnum i))
+ (lambda (x)
+ (when (< i glen)
+ (rplaca x (aref cseq (aref act i)))
+ (incf i))))
+ seq)))))))
+
+(defun permute (seq perm)
+ (declare (type sequence seq)
+ (type permutation perm))
+ (let ((cseq (copy-seq seq)))
+ (permute! cseq perm)))
+;;
+#+nil
+(defun permute-argument (func-symbol perm)
+ (declare (type symbol func-symbol)
+ (type permutation perm))
+ (let* ((glen (group-rank perm))
+ (args (loop for i from 0 below glen
+ collect (gensym))))
+ (eval `(lambda (,@args &rest rest)
+ (apply ',func-symbol (append (list ,@(permute! args perm)) rest))))))
+
+(defun argument-permute (func perm)
+ (declare (type function func)
+ (type permutation perm))
+ (lambda (&rest args)
+ (apply func (permute! args perm))))
+
+(defun curry (func perm &rest curried-args)
+ (declare (type function func)
+ (type permutation perm))
+ (lambda (&rest args)
+ (apply func (permute! (append curried-args args) perm))))
+
+(defun compose (func-a func-b perm)
+ (declare (type function func-a func-b)
+ (type permutation perm))
+ (lambda (&rest args)
+ (apply func-a (permute! (multiple-value-list (funcall func-b args)) perm))))
+;;
+
+(defun idx-sort-permute (seq predicate)
+ "
+ (sort-permute seq predicate)
+
+ Sorts a index-array and also returns
+ the permutation-action required to move
+ from the given sequence to the sorted form.
+
+ Takes about 10x the running time which can be
+ achieved with cl:sort.
+ "
+ (declare (type (index-array *) seq)
+ (type function predicate))
+ (let* ((len (length seq))
+ (perm (id-action-repr len)))
+ (declare (type index-type len)
+ (type (index-array *) perm))
+ (labels ((qsort-bounds (lb ub)
+ (declare (type index-type lb ub))
+ #+nil(format t "~a lb:~a ub:~a ~%" seq lb ub)
+ (if (= ub (1+ lb)) t
+ (let* ((ele (aref seq lb))
+ (ele-idx (very-quickly
+ (loop
+ for i of-type index-type from (1+ lb) below ub
+ with ele-idx of-type index-type = lb
+ do (unless (funcall predicate ele (aref seq i))
+ (when (> i (1+ ele-idx))
+ (rotatef (aref seq ele-idx) (aref seq (1+ ele-idx)))
+ (rotatef (aref perm ele-idx) (aref perm (1+ ele-idx))))
+ (rotatef (aref seq ele-idx) (aref seq i))
+ (rotatef (aref perm ele-idx) (aref perm i))
+ (incf ele-idx)
+ #+nil(format t " ~a ~%" seq))
+ finally (return ele-idx)))))
+ (when (> (- ub ele-idx) 2)
+ (qsort-bounds (1+ ele-idx) ub))
+ (when (> (- ele-idx lb) 1)
+ (qsort-bounds lb ele-idx))))))
+ (qsort-bounds 0 len)
+ (values seq (action->cycle (make-paction perm))))))
diff --git a/src/permutations.lisp b/src/permutations.lisp
deleted file mode 100644
index c68d97d..0000000
--- a/src/permutations.lisp
+++ /dev/null
@@ -1,230 +0,0 @@
-(in-package :matlisp)
-
-(define-condition permutation-error (generic-error)
- ((message :reader message :initform "Object is not a permutation."))
- (:documentation "Object is not a permutation."))
-
-;;Class definitions----------------------------------------------;;
-(defclass permutation ()
- ((representation :accessor repr
- :initarg :repr)
- (group-rank :accessor group-rank
- :type index-type)))
-;;
-(defclass permutation-cycle (permutation)
- ((representation :type cons)))
-
-(defmethod initialize-instance :after ((per permutation-cycle) &rest initargs)
- (declare (ignore initargs))
- (let ((cls 0))
- (declare (type index-type cls))
- (unless (very-quickly
- (dolist (cyc (r-value per) t)
- (unless (cycle-p cyc)
- (return nil))
- (setf cls (max cls (idx-max cyc)))))
- (error 'permutation-error))
- (setf (group-rank per) (the index-type (1+ cls)))))
-;;
-(defclass permutation-action (permutation)
- ((representation :type (index-array *))))
-
-(defmethod initialize-instance :after ((per permutation-action) &rest initargs)
- (declare (ignore initargs))
- (let ((act (r-value per)))
- (declare (type (index-array *) act))
- (unless (action-p act)
- (error 'permutation-error))
- (setf (group-rank per) (idx-max act))))
-
-;;Conversions and validation-------------------------------------;;
-(defun insert-element (x sort l-b u-b)
- "Does a binary-esque sort to keep track of elements in
- a permutation, in descending order. If there are duplicates
- of X in sort between L-B and U-B (both inclusive), or if X < 0,
- then throws a PERMUTATION-ERROR."
- (declare (type index-type x l-b u-b)
- (type (index-array *) sort))
- (let* ((len u-b))
- (labels ((insert-ele (l-b u-b)
- (declare (type index-type l-b u-b))
- (let* ((midx (+ l-b (floor (- u-b l-b) 2)))
- (mid (aref sort midx)))
- (declare (type index-type midx mid))
- (cond
- ((or (< x 0) (member x `(,(aref sort u-b) ,(aref sort l-b) ,mid)))
- (error 'permutation-error))
- ((= midx l-b)
- (when (> x (aref sort u-b))
- (very-quickly
- (loop
- with sidx of-type index-type = (+ midx (if (> x mid) 0 1))
- for i of-type index-type downfrom (1- len) to sidx
- do (setf (aref sort (+ i 1)) (aref sort i))
- finally (setf (aref sort sidx) x)))))
- ((< x mid) (insert-ele midx u-b))
- ((> x mid) (insert-ele l-b midx)))
- sort)))
- (insert-ele l-b u-b))))
-
-(defun cycle-new-p (perm)
- "Does a sorting operation to check for duplicate elements in
- the cycle representation of a permutation."
- (declare (type (index-array *) perm))
- (let* ((len (length perm))
- (sort (very-quickly (sort (copy-seq perm) #'<))))
- (declare (type (index-array *) sort)
- (type index-type len))
- (very-quickly
- (loop for i of-type index-type from 1 below len
- when (= (aref sort i) (aref sort (1- i)))
- do (return nil)
- finally (return t)))))
-
-(defun action-p (act)
- "Checks if ARR is a possible permutation vector. A permutation pi
- is characterized by a vector containing the indices from 0,...,
- @function{length}(@arg{perm})-1 in some order."
- (declare (type (index-array *) act))
- (let* ((len (length act))
- (sort (very-quickly (sort (copy-seq act) #'<))))
- (declare (type (index-array *) sort)
- (type index-type len))
- (very-quickly
- (loop for i of-type index-type from 0 below len
- unless (= (aref sort i) i)
- do (return nil)
- finally (return t)))))
-
-(defun action->cycle (act)
- ;;Caution: will go into an infinite loop if object is not proper.
- "
- This function obtains the canonical cycle representation
- of a permutation. The first argument is the action of the
- permutation on the array #(0 1 2 3 ..).
- \"Canonical\" may be a bit of an overstatement; this is the way
- S_n was presented by Van der Waerden.
-"
- (declare (type permutation-action per))
- (mlet*
- ((arr (r-value per) :type (index-array *)))
- (labels ((find-cycle (arr x0)
- "This function obtains a permutation cycle starting from x_0.
- The first argument is the action of the permutation on the
- array #(0 1 2 ..)"
- (declare (type (index-array *) arr)
- (type index-type x0))
- (if (= (aref arr x0) x0) (values #() nil)
- (destructuring-bind (n lst)
- (do ((i 0 (+ i 1))
- (x x0 (aref arr x))
- (ret nil (cons x ret))
- (count 0 (+ count (if (= x x0) 1 0))))
- ((and (= count 1) (= x x0)) (list i ret)))
- (values (make-array n :element-type 'index-type :initial-contents lst) lst))))
- (cycle-walk (cyc ignore)
- (declare (optimize (speed 3) (safety 0)))
- (let ((x0 (find-if-not #'(lambda (x) (member x ignore)) arr)))
- (if (null x0) cyc
- (multiple-value-bind (cnew clst) (find-cycle arr x0)
- (cycle-walk (if (null clst) cyc (cons cnew cyc))
- (nconc ignore (if (null clst) (list x0) clst))))))))
- (cycle-walk nil nil))))
-;;---------------------------------------------------------------;;
-
-
-(defun cycles->action (cyc)
- )
-
-;;
-(defun apply-cycle! (seq cyc)
- (declare (type (index-array *) cyc)
- (type (vector * *) seq))
- (unless (cycle-p cyc)
- (error 'permutation-error))
- (when (> (length cyc) 1)
- (let ((xl (aref seq (aref cyc (- (length cyc) 1)))))
- (loop for i downfrom (- (length cyc) 1) to 0
- do (setf (aref seq (aref cyc i))
- (if (= i 0) xl
- (aref seq (aref cyc (- i 1))))))))
- seq)
-
-(defun permute! (seq cycs)
- (unless (or (null cycs) (= (length seq) 0))
- (dolist (cyc cycs)
- (apply-cycle! seq cyc)))
- seq)
-
-(defun arg-perm (func cycs)
- (if (null cycs)
- func
- (lambda (&rest args)
- (let ((argvec (permute! (apply #'vector args) cycs)))
- (apply func (loop for i from 0 below (length argvec)
- collect (aref argvec i)))))))
-
-(defun compose (func func)
-
-;; (defun compose (..)
-;; )
-
-(defun seqrnd (seq)
- "Randomize the elements of a sequence. Destructive on SEQ."
- (sort seq #'> :key #'(lambda (x) (random 1.0))))
-
-;;
-
-(defun allocate-unit-permutation (n)
- (declare (type fixnum n))
- (let ((ret (allocate-index-store n)))
- (declare (type (index-array *) ret))
- (very-quickly
- (loop
- for i of-type index-type from 0 below n
- do (setf (aref ret i) i)))
- ret))
-
-(defun sort-permute (seq predicate)
- "
- (sort-permute seq predicate)
-
- Sorts a index-array and also returns
- the permutation-action required to move
- from the given sequence to the sorted form.
-
- Takes about 10x the running time which can be
- achieved with cl:sort.
- "
- (declare (type (index-array *) seq)
- (type function predicate))
- (let* ((len (length seq))
- (perm (allocate-unit-permutation len)))
- (declare (type index-type len)
- (type (index-array *) perm))
- (labels ((qsort-bounds (lb ub)
- (declare (type index-type lb ub))
- #+nil(format t "~a lb:~a ub:~a ~%" seq lb ub)
- (if (= ub (1+ lb)) t
- (let* ((ele (aref seq lb))
- (ele-idx (very-quickly
- (loop
- for i of-type index-type from (1+ lb) below ub
- with ele-idx of-type index-type = lb
- do (unless (funcall predicate ele (aref seq i))
- (when (> i (1+ ele-idx))
- (rotatef (aref seq ele-idx) (aref seq (1+ ele-idx)))
- (rotatef (aref perm ele-idx) (aref perm (1+ ele-idx))))
- (...
[truncated message content] |
Thanks for helping keep SourceForge clean.
X