From: Akshay S. <ak...@us...> - 2013-02-02 18:30:51
|
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 7a8ab5c0db938424bfc8d2ebca6022e2673b6a9a (commit) via 28177ada352e525b622e38a5c7baf97986854288 (commit) from 702e41559e7f9cc9d84d108e467a1ec132e6f50f (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 7a8ab5c0db938424bfc8d2ebca6022e2673b6a9a Author: Akshay Srinivasan <aks...@gm...> Date: Sat Feb 2 10:25:32 2013 -0800 Saving changes in permutation.lisp diff --git a/packages.lisp b/packages.lisp index 6e70c2a..87b3afc 100644 --- a/packages.lisp +++ b/packages.lisp @@ -66,7 +66,7 @@ (defpackage "MATLISP-UTILITIES" (:use #:common-lisp #:matlisp-conditions) - (:export #:ensure-list + (:export #:ensure-list #:id ;;#:make-ring #:zip #:zip-eq #:cut-cons-chain! #:slot-values diff --git a/src/base/permutation.lisp b/src/base/permutation.lisp index 9439689..9f5bbda 100644 --- a/src/base/permutation.lisp +++ b/src/base/permutation.lisp @@ -1,35 +1,33 @@ (in-package #:matlisp) ;;This must match the type used in LAPACK -(deftype perrepr-type () +(deftype pindex-type () '(unsigned-byte 32)) -(deftype perrepr-vector (&optional (size '*)) - `(simple-array perrepr-type (,size))) +(deftype pindex-store-vector (&optional (size '*)) + `(simple-array pindex-type (,size))) -(make-array-allocator allocate-perrepr-store 'perrepr-type 0 - " +(make-array-allocator allocate-pindex-store 'pindex-type 0 + " Syntax ====== - (ALLOCATE-PERREPR-STORE SIZE [INITIAL-ELEMENT 0]) + (ALLOCATE-PINDEX-STORE SIZE [INITIAL-ELEMENT 0]) Purpose ======= Allocates integer4 (32-bits) storage.") ;; -(defun perrepr-id-action (n) +(definline pindex-id (n) (declare (type fixnum n)) - (let ((ret (allocate-perrepr-store n))) - (declare (type perrepr-vector ret)) - (very-quickly - (loop - :for i :of-type perrepr-type :from 0 :below n - :do (setf (aref ret i) i))) - ret)) + (let-typed ((ret (allocate-pindex-store n) :type pindex-store-vector)) + (very-quickly + (loop :for i :of-type pindex-type :from 0 :below n + :do (setf (aref ret i) i))) + ret)) -(definline perv (&rest contents) - (make-array (length contents) :element-type 'perrepr-type :initial-contents contents)) +(definline pidxv (&rest contents) + (make-array (length contents) :element-type 'pindex-type :initial-contents contents)) ;;Write a uniform randomiser (defun seqrnd (seq) @@ -39,141 +37,77 @@ ;;Class definitions----------------------------------------------;; (defclass permutation () - ((representation :accessor repr - :initarg :repr) - (group-rank :accessor group-rank - :type index-type))) + ((store :accessor store :initarg :store) + (permutation-size :accessor permutation-size :type index-type))) (defmethod print-object ((per permutation) stream) (print-unreadable-object (per stream :type t) - (format stream "S_~a~%" (group-rank per)) - (if (zerop (group-rank per)) + (format stream "S_~a~%" (permutation-size per)) + (if (= (permutation-size per) 1) (format stream "ID~%") - (format stream "~a~%" (repr per))))) -;; + (format stream "~a~%" (store per))))) +;; (defclass permutation-action (permutation) - ((representation :type perrepr-vector))) + ((store :type pindex-store-vector))) - -(defmethod initialize-instance :after ((per permutation-action) &rest initargs) +(defmethod initialize-instance :after ((perm permutation-action) &rest initargs) (declare (ignore initargs)) - (labels ((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 'perrepr-vector)) nil - (locally - (declare (type perrepr-vector act)) - (let* ((len (length act)) - (sort (very-quickly (sort (copy-seq act) #'<)))) - (declare (type perrepr-vector 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)))))))) - (assert (action-repr-p (repr per)) nil 'permutation-invalid-error) - (setf (group-rank per) (length (repr per))))) -;; + (when *check-after-initializing?* + (let-typed ((repr (store perm) :type pindex-store-vector)) + (very-quickly + (loop :for i :of-type index-type :from 0 :below (length repr) + :with srepr :of-type pindex-store-vector := (sort (copy-seq repr) #'<) + :do (assert (= (aref srepr i) i) nil 'permutation-invalid-error))) + (setf (permutation-size perm) (length repr))))) +;; (defclass permutation-cycle (permutation) - ((representation :type list))) + ((store :type list))) (defmethod initialize-instance :after ((per permutation-cycle) &rest initargs) (declare (ignore initargs)) - (labels ((cycle-repr-p (perm) - " - Does a sorting operation to check for duplicate elements in - the cycle representation of a permutation. -" - (if (not (typep perm 'perrepr-vector)) nil - (locally - (declare (type perrepr-vector perm)) - (let ((len (length perm))) - (declare (type index-type len)) - (if (<= len 1) nil - (let ((sort (very-quickly (sort (copy-seq perm) #'<)))) - (declare (type perrepr-vector 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)))))))))) - (if (null (repr per)) (setf (group-rank per) 0) + (when *check-after-initializing?* + (if (null (store per)) + (setf (permutation-size per) 1) (very-quickly (loop - :for cyc :of-type perrepr-vector :in (repr per) - :unless (cycle-repr-p cyc) - :do (error 'permutation-invalid-error) - :maximizing (lvec-max cyc) into g-rnk of-type index-type - :finally (setf (group-rank per) (the index-type (1+ g-rnk)))))))) + :for cyc :of-type pindex-store-vector :in (store per) + :with ss :of-type pindex-type := 0 + :do (loop + :for i :of-type index-type :from 1 :below (length cyc) + :with scyc :of-type pindex-store-vector := (sort (copy-seq cyc) #'<) + :do (assert (/= (aref scyc (1- i)) (aref scyc i)) nil 'permutation-invalid-error) + :finally (setf ss (max ss (aref scyc (1- (length scyc)))))) + :finally (setf (permutation-size per) (1+ ss))))))) ;; (defclass permutation-pivot-flip (permutation) - ((representation :type perrepr-vector))) + ((store :type pindex-store-vector))) (defmethod initialize-instance :after ((per permutation-pivot-flip) &rest initargs) (declare (ignore initargs)) - (labels ((pivot-flip-p (idiv) - " - Checks if ARR is a possible permutation pivot-flip. A pivot - flip is more algorithmic in its representation. If a sequence - is given, apply a pivot-flip on it is equivalent to running from - the left to the right of the permutation by flipping (pi(i), i) - sequentially. This is the representation used in LAPACK. -" - (if (not (typep idiv 'perrepr-vector)) nil - (let ((len (length idiv))) - (declare (type perrepr-vector idiv) - (type index-type len)) - (very-quickly - (loop :for i :of-type index-type :from 0 :below len - :do (unless (< -1 (aref idiv i) len) - (return nil)) - :finally (return t))))))) - (assert (pivot-flip-p (repr per)) nil 'permutation-invalid-error) - (setf (group-rank per) (length (repr per))))) - -;; -;; (defclass permutation-matrix (permutation) -;; ((representation :type perrepr-vector))) - -;; (defmethod initialize-instance :after ((per permutation-pivot-flip) &rest initargs) -;; (declare (ignore initargs)) -;; (labels ((pivot-flip-p (idiv) -;; " -;; Checks if ARR is a possible permutation pivot-flip. A pivot -;; flip is more algorithmic in its representation. If a sequence -;; is given, apply a pivot-flip on it is equivalent to running from -;; the left to the right of the permutation by flipping (pi(i), i) -;; sequentially. This is the representation used in LAPACK. -;; " -;; (if (not (typep idiv 'perrepr-vector)) nil -;; (let ((len (length idiv))) -;; (declare (type perrepr-vector idiv) -;; (type index-type len)) -;; (very-quickly -;; (loop :for i :of-type index-type :from 0 :below len -;; :do (unless (< -1 (aref idiv i) len) -;; (return nil)) -;; :finally (return t))))))) -;; (assert (pivot-flip-p (repr per)) nil 'permutation-invalid-error) -;; (setf (group-rank per) (length (repr per))))) - + (when *check-after-initializing?* + (let*-typed ((repr (store per) :type pindex-store-vector) + (len (length repr) :type index-type)) + (very-quickly + (loop :for i :of-type index-type :from 0 :below len + :do (assert (< -1 (aref repr i) len) nil 'permutation-invalid-error))) + (setf (store-size per) len)))) ;; -(definline make-pcycle (&rest args) - (make-instance 'permutation-cycle :repr args)) +(defclass permutation-matrix (permutation) + ((store :type pindex-store-vector))) -(definline make-paction (pact) - (make-instance 'permutation-action :repr pact)) - -(definline make-pidx (pact) - (make-instance 'permutation-pivot-flip :repr pact)) +(defmethod initialize-instance :after ((perm permutation-matrix) &rest initargs) + (declare (ignore initargs)) + (when *check-after-initializing?* + (let-typed ((repr (store perm) :type pindex-store-vector)) + (very-quickly + (loop :for i :of-type index-type :from 0 :below (length repr) + :with srepr :of-type pindex-store-vector := (sort (copy-seq repr) #'<) + :do (assert (= (aref srepr i) i) nil 'permutation-invalid-error))) + (setf (permutation-size perm) (length repr))))) ;;Generic permute! method. (defgeneric permute! (thing permutation &optional argument) @@ -195,311 +129,323 @@ (definline permute (thing perm &optional (arg 0)) (permute! (copy thing) perm arg)) -;;Action -(defmethod permute! ((seq cons) (perm permutation-action) &optional arg) - (declare (ignore arg)) - (let ((cseq (make-array (length seq) :initial-contents seq)) - (act (repr perm)) - (glen (group-rank perm))) - (mapl - (let ((i 0)) - (declare (type fixnum i)) - (lambda (x) - (when (< i glen) - (rplaca x (aref cseq (aref act i))) - (incf i)))) seq))) - -(defmethod permute! ((seq vector) (perm permutation-action) &optional arg) - (declare (ignore arg)) - (let ((cseq (make-array (length seq) :initial-contents seq)) - (act (repr perm))) - (loop - :for i :from 0 :below (group-rank perm) - :do (unless (= i (aref act i)) - (setf (aref seq i) (aref cseq (aref act i)))) - :finally (return seq)))) - -(defmethod permute! ((ten standard-tensor) (perm permutation-action) &optional (arg 0)) - (let ((cyc (action->cycle perm))) - (permute! ten cyc arg))) - -;;Cycle -;;Might be useful ? -(defun apply-cycle! (seq pcyc) - (declare (type perrepr-vector pcyc) - (type vector seq)) - (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)))))))) - -(defmethod permute! ((seq cons) (perm permutation-cycle) &optional arg) - (declare (ignore arg)) - (let ((cseq (make-array (length seq) :initial-contents seq)) - (glen (group-rank perm))) - (dolist (cyc (repr perm)) - (declare (type perrepr-vector cyc)) - (apply-cycle! cseq cyc)) - (mapl - (let ((i 0)) - (lambda (x) - (when (< i glen) - (rplaca x (aref cseq i)) - (incf i)))) seq))) - -(defmethod permute! ((seq vector) (perm permutation-cycle) &optional arg) - (declare (ignore arg)) - (dolist (cyc (repr perm) seq) - (declare (type perrepr-vector cyc)) - (apply-cycle! seq cyc))) - -(defmethod permute! ((A standard-tensor) (perm permutation-cycle) &optional (arg 0)) - (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) - (rplaca (nthcdr arg slst) 0) - (values (sub-tensor~ A slst) (sub-tensor~ A slst))) - (let-typed ((cyclst (repr perm) :type cons) - (cp-ten (make-instance (class-of tone) - :dimensions (copy-seq (dimensions tone)))) - (std-arg (aref (strides A) arg) :type index-type) - (hd-sl (head ttwo) :type index-type)) - (dolist (cyc cyclst) - (declare (type perrepr-vector cyc)) - (setf (head tone) (+ hd-sl (* std-arg (aref cyc (1- (length cyc)))))) - (copy! tone cp-ten) - (loop :for i :of-type index-type :downfrom (1- (length cyc)) :to 0 - :do (progn - (setf (head tone) (+ hd-sl (* std-arg (aref cyc i)))) - (copy! - (if (= i 0) cp-ten - (progn - (setf (head ttwo) (+ hd-sl (* std-arg (aref cyc (1- i))))) - ttwo)) - tone)))))) - A) - -;;Pivot idx -(defmethod permute! ((seq vector) (perm permutation-pivot-flip) &optional arg) - (declare (ignore arg)) - (let-typed ((pidx (repr perm) :type perrepr-vector)) - (loop :for i :of-type index-type :from 0 :below (group-rank perm) - :unless (= i (aref pidx i)) - :do (rotatef (aref seq i) (aref seq (aref pidx i))) - :finally (return seq)))) - -(defmethod permute! ((seq cons) (perm permutation-pivot-flip) &optional arg) - (declare (ignore arg)) - (let ((cseq (make-array (length seq) :initial-contents seq)) - (glen (group-rank perm))) - (permute! cseq perm) - (mapl - (let ((i 0)) - (lambda (x) - (when (< i glen) - (rplaca x (aref cseq i)) - (incf i)))) seq))) - -(defmethod permute! ((A standard-tensor) (perm permutation-pivot-flip) &optional (arg 0)) - (let ((idiv (repr perm))) - (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) - (rplaca (nthcdr arg slst) 0) - (values (sub-tensor~ A slst nil) (sub-tensor~ A slst nil))) - (let ((argstd (aref (strides A) arg)) - (hd-sl (head ttwo))) - (loop :for i :from 0 :below (length idiv) - :do (progn - (unless (= i (aref idiv i)) - (setf (head ttwo) (+ hd-sl (* (aref idiv i) argstd))) - (swap! tone ttwo)) - (incf (head tone) argstd)))))) - A) - -;;Conversions----------------------------------------------------;; -(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 perrepr-vector)) - (labels ((find-cycle (x0) - ;; This function obtains the cycle starting from x_0. - (declare (type perrepr-type x0)) - (if (= (aref arr x0) x0) (values 0 nil) - (very-quickly - (loop - :for x :of-type perrepr-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 'perrepr-type :initial-contents clst) cyc)) - (nconc ignore (if (= clen 0) (list x0) clst)))))))) - (let ((cyc-lst (cycle-walk nil nil))) - (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 (perrepr-id-action (group-rank cyc))) - (cycs-repr (repr cyc))) - (declare (type perrepr-vector act-repr)) - (dolist (cyc cycs-repr) - (declare (type perrepr-vector 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))) - -(defun pivot-flip->action (pflip) - (declare (type permutation-pivot-flip pflip)) - (let* ((idiv (repr pflip)) - (len (length idiv))) - (declare (type perrepr-vector idiv) - (type index-type len)) - (let ((act (perrepr-id-action len))) - (declare (type perrepr-vector act)) - (very-quickly - (loop :for i :from 0 :below len - :do (let ((val (aref idiv i))) - (unless (= val i) - (rotatef (aref act i) (aref act val)))))) - (make-instance 'permutation-action :repr act)))) - -(defun mod-max (seq lidx uidx) - (declare (type perrepr-vector seq)) - (let ((len (length seq))) - (very-quickly - (loop :for idx :of-type index-type :downfrom uidx :above lidx - :with max :of-type perrepr-type := (aref seq uidx) - :with max-idx :of-type index-type := uidx - :do (let ((ele (aref seq (mod idx len)))) - (when (> ele max) - (setf max ele - max-idx idx))) - :finally (return (values max max-idx)))))) - -#| - -(defun cycle->pivot-flip (cyc) - (let ((cp-cyc (copy-seq cyc))) - (let - (labels ((mod-max (seq lidx uidx) - (declare (type perrepr-vector seq)) - (let ((len (length cyc))) - (very-quickly - (loop :for idx :of-type index-type :downfrom uidx :above lidx - :with max :of-type perrepr-type := (aref seq uidx) - :with max-idx :of-type index-type := uidx - :do (let ((ele (aref seq (mod idx len)))) - (when (> ele max) - (setf max ele - max-idx idx))) - :finally (return (values max max-idx)))))) - (get-flip (lidx uidx) - (multiple-value-bind (max max-idx) (mod-max cyc lidx uidx) - (let ((ele-0 (aref cyc (mod max-idx len))) - (ele-1 (aref cyc (mod (1- max-idx) len)))) - (setf (aref pidx (max ele-0 ele-1)) - (min ele-0 ele-1)) - |# - -#+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)))) -;; - -(definline sort-permute (seq predicate) +;; ;;Action +;; (defmethod permute! ((seq cons) (perm permutation-action) &optional arg) +;; (declare (ignore arg)) +;; (let ((cseq (make-array (length seq) :initial-contents seq)) +;; (act (repr perm)) +;; (glen (group-rank perm))) +;; (mapl +;; (let ((i 0)) +;; (declare (type fixnum i)) +;; (lambda (x) +;; (when (< i glen) +;; (rplaca x (aref cseq (aref act i))) +;; (incf i)))) seq))) + +;; (defmethod permute! ((seq vector) (perm permutation-action) &optional arg) +;; (declare (ignore arg)) +;; (let ((cseq (make-array (length seq) :initial-contents seq)) +;; (act (repr perm))) +;; (loop +;; :for i :from 0 :below (group-rank perm) +;; :do (unless (= i (aref act i)) +;; (setf (aref seq i) (aref cseq (aref act i)))) +;; :finally (return seq)))) + +;; (defmethod permute! ((ten standard-tensor) (perm permutation-action) &optional (arg 0)) +;; (let ((cyc (action->cycle perm))) +;; (permute! ten cyc arg))) + +;; ;;Cycle +;; ;;Might be useful ? +;; (defun apply-cycle! (seq pcyc) +;; (declare (type pindex-store-vector pcyc) +;; (type vector seq)) +;; (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)))))))) + +;; (defmethod permute! ((seq cons) (perm permutation-cycle) &optional arg) +;; (declare (ignore arg)) +;; (let ((cseq (make-array (length seq) :initial-contents seq)) +;; (glen (group-rank perm))) +;; (dolist (cyc (repr perm)) +;; (declare (type pindex-store-vector cyc)) +;; (apply-cycle! cseq cyc)) +;; (mapl +;; (let ((i 0)) +;; (lambda (x) +;; (when (< i glen) +;; (rplaca x (aref cseq i)) +;; (incf i)))) seq))) + +;; (defmethod permute! ((seq vector) (perm permutation-cycle) &optional arg) +;; (declare (ignore arg)) +;; (dolist (cyc (repr perm) seq) +;; (declare (type pindex-store-vector cyc)) +;; (apply-cycle! seq cyc))) + +;; (defmethod permute! ((A standard-tensor) (perm permutation-cycle) &optional (arg 0)) +;; (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) +;; (rplaca (nthcdr arg slst) 0) +;; (values (sub-tensor~ A slst) (sub-tensor~ A slst))) +;; (let-typed ((cyclst (repr perm) :type cons) +;; (cp-ten (make-instance (class-of tone) +;; :dimensions (copy-seq (dimensions tone)))) +;; (std-arg (aref (strides A) arg) :type index-type) +;; (hd-sl (head ttwo) :type index-type)) +;; (dolist (cyc cyclst) +;; (declare (type pindex-store-vector cyc)) +;; (setf (head tone) (+ hd-sl (* std-arg (aref cyc (1- (length cyc)))))) +;; (copy! tone cp-ten) +;; (loop :for i :of-type index-type :downfrom (1- (length cyc)) :to 0 +;; :do (progn +;; (setf (head tone) (+ hd-sl (* std-arg (aref cyc i)))) +;; (copy! +;; (if (= i 0) cp-ten +;; (progn +;; (setf (head ttwo) (+ hd-sl (* std-arg (aref cyc (1- i))))) +;; ttwo)) +;; tone)))))) +;; A) + +;; ;;Pivot idx +;; (defmethod permute! ((seq vector) (perm permutation-pivot-flip) &optional arg) +;; (declare (ignore arg)) +;; (let-typed ((pidx (repr perm) :type pindex-store-vector)) +;; (loop :for i :of-type index-type :from 0 :below (group-rank perm) +;; :unless (= i (aref pidx i)) +;; :do (rotatef (aref seq i) (aref seq (aref pidx i))) +;; :finally (return seq)))) + +;; (defmethod permute! ((seq cons) (perm permutation-pivot-flip) &optional arg) +;; (declare (ignore arg)) +;; (let ((cseq (make-array (length seq) :initial-contents seq)) +;; (glen (group-rank perm))) +;; (permute! cseq perm) +;; (mapl +;; (let ((i 0)) +;; (lambda (x) +;; (when (< i glen) +;; (rplaca x (aref cseq i)) +;; (incf i)))) seq))) + +;; (defmethod permute! ((A standard-tensor) (perm permutation-pivot-flip) &optional (arg 0)) +;; (let ((idiv (repr perm))) +;; (multiple-value-bind (tone ttwo) (let ((slst (make-list (rank A) :initial-element '\:))) +;; (rplaca (nthcdr arg slst) 0) +;; (values (sub-tensor~ A slst nil) (sub-tensor~ A slst nil))) +;; (let ((argstd (aref (strides A) arg)) +;; (hd-sl (head ttwo))) +;; (loop :for i :from 0 :below (length idiv) +;; :do (progn +;; (unless (= i (aref idiv i)) +;; (setf (head ttwo) (+ hd-sl (* (aref idiv i) argstd))) +;; (swap! tone ttwo)) +;; (incf (head tone) argstd)))))) +;; A) + +;; ;;Conversions----------------------------------------------------;; +;; (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 pindex-store-vector)) +;; (labels ((find-cycle (x0) +;; ;; This function obtains the cycle starting from x_0. +;; (declare (type pindex-type x0)) +;; (if (= (aref arr x0) x0) (values 0 nil) +;; (very-quickly +;; (loop +;; :for x :of-type pindex-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 'pindex-type :initial-contents clst) cyc)) +;; (nconc ignore (if (= clen 0) (list x0) clst)))))))) +;; (let ((cyc-lst (cycle-walk nil nil))) +;; (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 (pindex-id-action (group-rank cyc))) +;; (cycs-repr (repr cyc))) +;; (declare (type pindex-store-vector act-repr)) +;; (dolist (cyc cycs-repr) +;; (declare (type pindex-store-vector 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))) + +;; (defun pivot-flip->action (pflip) +;; (declare (type permutation-pivot-flip pflip)) +;; (let* ((idiv (repr pflip)) +;; (len (length idiv))) +;; (declare (type pindex-store-vector idiv) +;; (type index-type len)) +;; (let ((act (pindex-id-action len))) +;; (declare (type pindex-store-vector act)) +;; (very-quickly +;; (loop :for i :from 0 :below len +;; :do (let ((val (aref idiv i))) +;; (unless (= val i) +;; (rotatef (aref act i) (aref act val)))))) +;; (make-instance 'permutation-action :repr act)))) + +;; (defun mod-max (seq lidx uidx) +;; (declare (type pindex-store-vector seq)) +;; (let ((len (length seq))) +;; (very-quickly +;; (loop :for idx :of-type index-type :downfrom uidx :above lidx +;; :with max :of-type pindex-type := (aref seq uidx) +;; :with max-idx :of-type index-type := uidx +;; :do (let ((ele (aref seq (mod idx len)))) +;; (when (> ele max) +;; (setf max ele +;; max-idx idx))) +;; :finally (return (values max max-idx)))))) + +;; #| + +;; (defun cycle->pivot-flip (cyc) +;; (let ((cp-cyc (copy-seq cyc))) +;; (let +;; (labels ((mod-max (seq lidx uidx) +;; (declare (type pindex-store-vector seq)) +;; (let ((len (length cyc))) +;; (very-quickly +;; (loop :for idx :of-type index-type :downfrom uidx :above lidx +;; :with max :of-type pindex-type := (aref seq uidx) +;; :with max-idx :of-type index-type := uidx +;; :do (let ((ele (aref seq (mod idx len)))) +;; (when (> ele max) +;; (setf max ele +;; max-idx idx))) +;; :finally (return (values max max-idx)))))) +;; (get-flip (lidx uidx) +;; (multiple-value-bind (max max-idx) (mod-max cyc lidx uidx) +;; (let ((ele-0 (aref cyc (mod max-idx len))) +;; (ele-1 (aref cyc (mod (1- max-idx) len)))) +;; (setf (aref pidx (max ele-0 ele-1)) +;; (min ele-0 ele-1)) +;; |# + +;; #+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)))) +;; ;; + + +;; (defstruct (ring (:constructor nil)) +;; (circular-list nil :type list) +;; (end-list nil :type list)) + +;; (defun make-ring (n) +;; (let ((ret (cons nil nil))) +;; (loop :for i :from 0 :below (1- n) +;; :with tail :of-type cons := ret +;; :do (setf (cdr tail) (cons nil nil) +;; tail (cdr tail)) +;; :finally (setf (cdr tail) ret)) +;; ret)) + + +;;This function is ugly of-course, but is also very very quick! +(definline sort-permute (seq predicate &key (key #'matlisp-utilities:id)) " Sorts a lisp-vector in-place, by using the function @arg{predicate} as the order. Also computes the permutation which would sort the original sequence @arg{seq}. " (declare (type vector seq)) - ;;This function is ugly of-course, but is also very very quick! (let*-typed ((len (length seq) :type fixnum) - (perm (perrepr-id-action len) :type perrepr-vector) + (perm (pindex-id len) :type pindex-store-vector) (jobs (list `(0 ,len)))) - (loop ;;:repeat 10 - :for bounds := (pop jobs) :then (pop jobs) + (loop + :for bounds := (car jobs) :then (pop jobs) :until (null bounds) :do (let*-typed ((below-idx (first bounds) :type fixnum) (above-idx (second bounds) :type fixnum) (piv (+ below-idx (floor (- above-idx below-idx) 2)) :type fixnum)) - (loop ;;:repeat 10 - :with ele := (aref seq piv) + (loop + :with ele := (funcall key (aref seq piv)) :with lbound :of-type fixnum := below-idx :with ubound :of-type fixnum := (1- above-idx) :until (progn - ;;(format t "~%~a ~%" (list lbound piv ubound)) (loop :for i :of-type fixnum :from lbound :to piv - :until (or (= i piv) (funcall predicate ele (aref seq i))) + :until (or (= i piv) (funcall predicate ele (funcall key (aref seq i)))) :finally (setq lbound i)) (loop :for i :of-type fixnum :downfrom ubound :to piv - :until (or (= i piv) (funcall predicate (aref seq i) ele)) + :until (or (= i piv) (funcall predicate ele (funcall key (aref seq i)))) :finally (setq ubound i)) - ;;(format t "~a ~%" (list lbound piv ubound)) (cond ((= ubound lbound piv) (when (> (- piv below-idx) 1) (push `(,below-idx ,piv) jobs)) (when (> (- above-idx (1+ piv)) 1) (push `(,(1+ piv) ,above-idx) jobs)) - ;;(format t "~a~%" jobs) t) ((< lbound piv ubound) (rotatef (aref seq lbound) (aref seq ubound)) @@ -525,4 +471,4 @@ (decf piv) (decf ubound) nil))))) - :finally (return (values seq (make-instance 'permutation-action :repr perm)))))) + :finally (return (values seq (make-instance 'permutation-action :store perm)))))) diff --git a/src/utilities/functions.lisp b/src/utilities/functions.lisp index 079c6ba..3c5c40a 100644 --- a/src/utilities/functions.lisp +++ b/src/utilities/functions.lisp @@ -2,6 +2,19 @@ ;;These functions are used all over the place inside Matlisp's macros. (eval-when (:compile-toplevel :load-toplevel :execute) + + (declaim (inline id)) + (defun id (x) x) + + #+nil + (defun make-ring (n) + (let ((ret (cons nil nil))) + (loop :for i :from 0 :below (1- n) + :with tail :of-type cons := ret + :do (setf (cdr tail) (cons nil nil) + tail (cdr tail)) + :finally (setf (cdr tail) ret)) + ret)) (declaim (inline slot-values)) (defun slot-values (obj slots) diff --git a/src/utilities/macros.lisp b/src/utilities/macros.lisp index bad1160..259fb95 100644 --- a/src/utilities/macros.lisp +++ b/src/utilities/macros.lisp @@ -63,11 +63,12 @@ ,@body)))) (defmacro make-array-allocator (allocator-name type init &optional doc) - `(definline ,allocator-name (size &optional (initial-element ,init)) - ,@(unless (null doc) - `(,doc)) - (make-array size - :element-type ,type :initial-element initial-element))) + `(eval-when (:compile-toplevel :load-toplevel :execute) + (definline ,allocator-name (size &optional (initial-element ,init)) + ,@(unless (null doc) + `(,doc)) + (make-array size + :element-type ,type :initial-element initial-element)))) (defmacro let-typed (bindings &rest body) " commit 28177ada352e525b622e38a5c7baf97986854288 Author: Akshay Srinivasan <aks...@gm...> Date: Fri Feb 1 22:40:07 2013 -0800 Cosmetic changes to permutation.lisp, dlsode.lisp. diff --git a/src/base/permutation.lisp b/src/base/permutation.lisp index 0a69941..9439689 100644 --- a/src/base/permutation.lisp +++ b/src/base/permutation.lisp @@ -477,7 +477,6 @@ (loop ;;:repeat 10 :for bounds := (pop jobs) :then (pop jobs) :until (null bounds) - :finally (return (values seq (make-instance 'permutation-action :repr perm))) :do (let*-typed ((below-idx (first bounds) :type fixnum) (above-idx (second bounds) :type fixnum) (piv (+ below-idx (floor (- above-idx below-idx) 2)) :type fixnum)) @@ -525,4 +524,5 @@ (rotatef (aref perm lbound) (aref perm piv))) (decf piv) (decf ubound) - nil)))))))) + nil))))) + :finally (return (values seq (make-instance 'permutation-action :repr perm)))))) diff --git a/src/packages/odepack/dlsode.lisp b/src/packages/odepack/dlsode.lisp index a18b957..dab15bb 100644 --- a/src/packages/odepack/dlsode.lisp +++ b/src/packages/odepack/dlsode.lisp @@ -1,24 +1,11 @@ (in-package #:matlisp) (cffi:define-foreign-library libodepack - #+nil(:unix #.(translate-logical-pathname - (merge-pathnames "matlisp:lib;libodepack" - *shared-library-pathname-extension*))) (t (:default "libodepack"))) (cffi:use-foreign-library libodepack) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -#+nil -(def-fortran-routine testde :void - (field (:callback :void - (c-neq :integer :input) - (c-t :double-float :input) - (c-y (* :double-float :size c-neq) :input) - (c-ydot (* :double-float :size c-neq) :output))) - (neq :integer :input) - (y (* :double-float) :input-output)) - (def-fortran-routine dlsode :void "DLSODE in ODEPACK" (field (:callback dlsode-field-callback ----------------------------------------------------------------------- Summary of changes: packages.lisp | 2 +- src/base/permutation.lisp | 764 ++++++++++++++++++-------------------- src/packages/odepack/dlsode.lisp | 13 - src/utilities/functions.lisp | 13 + src/utilities/macros.lisp | 11 +- 5 files changed, 375 insertions(+), 428 deletions(-) hooks/post-receive -- matlisp |