Menu
▾
▴
matlisp-commit
|
From: Akshay S. <ak...@us...> - 2013-08-24 03:33:49
|
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, classy has been updated
via c6c440e0043ee6633cb729a0bc590e9ca97d5eff (commit)
via 6fb6102eb0a39b2bf48ba2dfbe98f2d7a1966935 (commit)
via 60f2b8490d5a31c90886c51d081bdb802f5431d6 (commit)
via 050548a939e93208a4990b6de248e3e39b3caa45 (commit)
via 7aa696273c37819ebc5b9ee1040d0f194dd8145a (commit)
via 570ae7eb80324580ee27cfa7ba1d20b11f779e41 (commit)
via 79a87b8605dba8ae97c8f354e42f1f081b127771 (commit)
via fc524fd099c95abfc3af0280e8a200e461cc9493 (commit)
via a3ea5898fd82c16b66bf0ce3d5615e370deb40a8 (commit)
via 69ca54a98c4e4a03e004268297644094b5541cae (commit)
via 6c6f96e88fab82f42a2cd563c53e90c48eb8da24 (commit)
via 1f3bc86e39b2b9d1a23946434486a99faf9d7eaf (commit)
via fd0546544cd3c21641688e03ff221b031ac01ae4 (commit)
via 62f126ef664982e0c8ffd132ade8bb5308833f56 (commit)
via 9afd000d1a6b497e3bf4fdc0318884b412773de7 (commit)
via a7bad8cec909a69bb312917406d3dfc1626f8c12 (commit)
via 831a7c79908907c3702b623d0eae3a0a1f746a58 (commit)
via 05cd941088d8c303f3b5f81d34a6fb336f1033a9 (commit)
via f051c33ed570af222cff1fbf93802dc8844034ad (commit)
via f314da645424005358b1921156218063d8ea64c2 (commit)
via f47066fa8877ca56d4d58f36a5a08515593b8a2f (commit)
via 1894d2f6b21c756ec5bb2ddd443e6d38a7087f2b (commit)
from 808353d428ddc07d365bf1de8abcc86f0179ee08 (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 c6c440e0043ee6633cb729a0bc590e9ca97d5eff
Author: Akshay Srinivasan <aks...@gm...>
Date: Fri Aug 23 20:33:50 2013 -0700
Added seq, map functions.
diff --git a/matlisp.asd b/matlisp.asd
index 44384f4..b4b9ac0 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -168,7 +168,9 @@
(:module "matlisp-special"
:pathname "special"
:depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3")
- :components ((:file "random")))
+ :components ((:file "random")
+ (:file "map")
+ (:file "seq")))
#+nil
(:module "matlisp-sugar"
:pathname "sugar"
diff --git a/src/base/generic-copy.lisp b/src/base/generic-copy.lisp
index 543fa5c..aa30346 100644
--- a/src/base/generic-copy.lisp
+++ b/src/base/generic-copy.lisp
@@ -78,6 +78,10 @@
y))
(copy! x y)))
+(defmethod copy! ((x cons) (y standard-tensor))
+ ;;You seriously weren't expecting efficiency were you :) ?
+ (let ((arr (make-array (list-dimensions x) :initial-contents x)))
+ (copy! arr y)))
;;
(defgeneric copy (object)
(:documentation
diff --git a/src/classes/numeric.lisp b/src/classes/numeric.lisp
index 6f3e87e..e266cfc 100644
--- a/src/classes/numeric.lisp
+++ b/src/classes/numeric.lisp
@@ -53,8 +53,9 @@
(deft/method t/l3-lb (sym complex-numeric-tensor) ()
'*complex-l3-fcall-lb*)
-;;SBCL uses specialized arrays for floating complex arrays.
-#-sbcl
+;;Comment this block if you want to use (simple-array (complex double-float) (*))
+;;as the underlying store. This will make Lisp-implementations of gemm .. faster
+;;but you'll lose the ability to use tensor-realpart~/imagpart~.
(progn
(deft/method t/store-element-type (sym complex-numeric-tensor) ()
(let ((cplx-type (macroexpand-1 `(t/field-type ,sym))))
diff --git a/src/level-1/map.lisp b/src/level-1/map.lisp
deleted file mode 100644
index 3359ea9..0000000
--- a/src/level-1/map.lisp
+++ /dev/null
@@ -1,72 +0,0 @@
-(in-package #:matlisp)
-
-(defgeneric mapsor! (func x y)
- (:documentation "
- y <- func(x)
-")
- (:method :before ((func function) (x standard-tensor) (y standard-tensor))
- (assert (very-quickly (lvec-eq (dimensions x) (dimensions y))) nil 'tensor-dimension-mismatch)))
-
-(defmethod mapsor! ((func function) (x standard-tensor) (y standard-tensor))
- (let ((clx (class-name (class-of x)))
- (cly (class-name (class-of y))))
- (assert (and
- (member clx *tensor-type-leaves*)
- (member cly *tensor-type-leaves*))
- nil 'tensor-abstract-class :tensor-class (list clx cly))
- (compile-and-eval
- `(defmethod mapsor! ((func function) (x ,clx) (y ,cly))
- (let ((sto-x (store x))
- (sto-y (store y)))
- (declare (type ,(store-type clx) sto-x)
- (type ,(store-type cly) sto-y))
- (very-quickly
- (mod-dotimes (idx (dimensions x))
- :with (linear-sums
- (of-x (strides x))
- (of-y (strides y)))
- :do (t/store-set ,cly (funcall func (t/store-ref ,clx sto-x of-x)) sto-y of-y))))
- y)))
- (mapsor! func x y))
-
-(defun mapcol (func x &optional (axis 1))
- (declare (type standard-tensor x))
- (assert (tensor-matrixp x) nil 'tensor-dimension-mismatch)
- (let* ((v-x (slice~ x axis))
- (st-x (aref (the index-store-vector (strides x)) axis))
- (ret (zeros (aref (the index-store-vector (dimensions x)) axis) (class-of x))))
- (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
- :do (progn
- (setf (ref ret i) (funcall func v-x))
- (incf (slot-value v-x 'head) st-x)))
- ret))
-
-(defun tensor-min (x)
- (let ((min nil))
- (mod-dotimes (idx (dimensions x))
- :do (when (or (null min) (< (ref x idx) min))
- (setf min (ref x idx))))
- min))
-
-(defun idx-min (v)
- (let ((min-idx nil)
- (min nil))
- (loop :for i :from 0 :below (aref (the index-store-vector (dimensions v)) 0)
- :do (when (or (null min) (< (ref v i) min))
- (setf min-idx i
- min (ref v i))))
- (values min-idx min)))
-
-(defun mapxis! (func x y &optional (axis 0))
- (declare (type standard-tensor x y))
- (multiple-value-bind (view-x view-y) (let ((slst (make-list (rank x) :initial-element '(* * *))))
- (rplaca (nthcdr axis slst) (list 0 '* 1))
- (values (sub-tensor~ x slst nil) (sub-tensor~ y slst nil)))
- (let ((st-x (aref (the index-store-vector (dimensions x)) axis))
- (st-y (aref (the index-store-vector (dimensions x)) axis)))
- (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
- :do (progn
- ;;May die if you're doing fancy stuff.
- (funcall func view-x view-y)
-
-
diff --git a/src/special/map.lisp b/src/special/map.lisp
new file mode 100644
index 0000000..9cb083a
--- /dev/null
+++ b/src/special/map.lisp
@@ -0,0 +1,70 @@
+(in-package #:matlisp)
+
+(defgeneric mapsor! (func x y)
+ (:documentation "
+ Syntax
+ ======
+ (MAPSOR! func x y)
+
+ Purpose
+ =======
+ Applies the function element-wise on x, and sets the corresponding
+ elements in y to the value returned by the function.
+
+ Example
+ =======
+ > (mapsor! #'sin (randn '(2 2)) (zeros '(2 2)))
+")
+ (:method :before ((func function) (x standard-tensor) (y standard-tensor))
+ (assert (very-quickly (lvec-eq (dimensions x) (dimensions y))) nil 'tensor-dimension-mismatch)))
+
+(defmethod mapsor! ((func function) (x standard-tensor) (y standard-tensor))
+ (let ((clx (class-name (class-of x)))
+ (cly (class-name (class-of y))))
+ (assert (and
+ (member clx *tensor-type-leaves*)
+ (member cly *tensor-type-leaves*))
+ nil 'tensor-abstract-class :tensor-class (list clx cly))
+ (compile-and-eval
+ `(defmethod mapsor! ((func function) (x ,clx) (y ,cly))
+ (let ((sto-x (store x))
+ (sto-y (store y)))
+ (declare (type ,(store-type clx) sto-x)
+ (type ,(store-type cly) sto-y))
+ (very-quickly
+ (mod-dotimes (idx (dimensions x))
+ :with (linear-sums
+ (of-x (strides x))
+ (of-y (strides y)))
+ :do (t/store-set ,cly (funcall func (t/store-ref ,clx sto-x of-x)) sto-y of-y))))
+ y)))
+ (mapsor! func x y))
+;;
+(defun mapslice (func x &optional (axis 1))
+ (declare (type standard-tensor x))
+ (let* ((v-x (slice~ x axis))
+ (st-x (aref (strides x) axis)))
+ (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
+ :collect (prog1 (funcall func (copy v-x))
+ (incf (slot-value v-x 'head) st-x)))))
+
+(defmacro tensor-foldl (type func ten init &optional (init-type (field-type type)))
+ (using-gensyms (decl (ten init))
+ (with-gensyms (sto idx of funcsym)
+ `(let* (,@decl
+ ,@(unless (symbolp func)
+ `((,funcsym ,func)))
+ (,sto (store ,ten)))
+ (declare (type ,type ,ten)
+ ,@(unless (symbolp func) `((type function ,funcsym)))
+ (type ,(store-type type) ,sto)
+ ,@(when init-type
+ `((type ,init-type ,init))))
+ (very-quickly
+ (mod-dotimes (,idx (dimensions ,ten))
+ :with (linear-sums
+ (,of (strides ,ten)))
+ :do (setf ,init (,@(if (symbolp func)
+ `(,func)
+ `(funcall ,funcsym)) ,init (t/store-ref ,type ,sto ,of)))))
+ ,init))))
diff --git a/src/sugar/seq.lisp b/src/special/seq.lisp
similarity index 70%
rename from src/sugar/seq.lisp
rename to src/special/seq.lisp
index 9c96efa..41ef589 100644
--- a/src/sugar/seq.lisp
+++ b/src/special/seq.lisp
@@ -27,22 +27,26 @@
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(in-package #:matlisp)
-(defun arange (start end &optional (h 1d0))
+(defun range (start end &optional (h 1d0))
(let ((quo (ceiling (if (> start end) (- start end) (- end start)) h)))
(if (= quo 0) nil
- (let*-typed ((ret (real-typed-zeros (idxv quo)) :type real-tensor)
- (sto-r (store ret) :type real-store-vector)
- (h (coerce-real-unforgiving (if (> start end) (- h) h)) :type real-type))
- (loop :for i :from 0 :below quo
- :for ori := (coerce-real-unforgiving start) :then (+ ori h)
- :do (setf (aref sto-r i) ori))
- ret))))
+ (let* ((ret (zeros quo 'real-tensor))
+ (sto (store ret))
+ (h (coerce (if (> start end) (- h) h) 'double-float)))
+ (declare (type (simple-array double-float (*)) sto)
+ (type double-float h))
+ (very-quickly
+ (loop :for i :from 0 :below quo
+ :for ori := (coerce start 'double-float) :then (+ ori h)
+ :do (t/store-set real-tensor ori sto i)))
+ ret))))
-(defun alinspace (start end &optional (num-points (1+ (abs (- start end)))))
- (let ((h (/ (- end start) (1- num-points))))
- (arange start (+ h end) (abs h))))
+(defun linspace (start end &optional (num-points (1+ (abs (- start end)))))
+ (let* ((num-points (floor num-points))
+ (h (/ (- end start) (1- num-points))))
+ (range start (+ h end) (abs h))))
-(defun range (start end &optional (h 1))
+(defun list-range (start end &optional (h 1))
(declare (type real start end h))
(let ((quo (ceiling (if (> start end) (- start end) (- end start)) h)))
(if (= quo 0) nil
@@ -51,6 +55,8 @@
:for ori := start :then (+ ori h)
:collect ori)))))
-(defun linspace (start end &optional (num-points (1+ (abs (- start end)))))
- (let ((h (/ (- end start) (1- num-points))))
- (range start (+ h end) (abs h))))
+(defun list-linspace (start end &optional (num-points (1+ (abs (- start end)))))
+ (let* ((num-points (floor num-points))
+ (h (/ (- end start) (1- num-points))))
+ (list-range start (+ h end) (abs h))))
+
diff --git a/src/utilities/lvec.lisp b/src/utilities/lvec.lisp
index 3ab1918..5e81a47 100644
--- a/src/utilities/lvec.lisp
+++ b/src/utilities/lvec.lisp
@@ -5,7 +5,7 @@
(declare (type vector))
(loop
:for i :of-type fixnum :from 0 :below (length vec)
- :for ret = (aref vec 0) :then (funcall func (aref vec i) ret)
+ :for ret = (aref vec 0) :then (funcall func ret (aref vec i))
:finally (return ret)))
(definline lvec-foldr (func vec)
commit 6fb6102eb0a39b2bf48ba2dfbe98f2d7a1966935
Author: Akshay Srinivasan <aks...@gm...>
Date: Thu Aug 22 21:58:42 2013 -0700
Cleaned up random.
diff --git a/AUTHORS b/AUTHORS
index 5f68edf..054d16e 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -9,3 +9,6 @@ Femlisp (www.femlisp.org); it has been used here
The infix reader is modified and included here with the
permission of its original author Mark Kantrowitz.
+
+The random number generators are borrowed and adapted from
+cl-random written by Tamas Papp.
\ No newline at end of file
diff --git a/matlisp.asd b/matlisp.asd
index cfb71d3..44384f4 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -165,6 +165,10 @@
:components ((:file "lu")
(:file "chol")
(:file "eig")))
+ (:module "matlisp-special"
+ :pathname "special"
+ :depends-on ("matlisp-base" "matlisp-classes" "matlisp-level-1" "matlisp-level-2" "matlisp-level-3")
+ :components ((:file "random")))
#+nil
(:module "matlisp-sugar"
:pathname "sugar"
diff --git a/src/packages/odepack/dlsode.lisp b/src/packages/odepack/dlsode.lisp
index dab15bb..466603c 100644
--- a/src/packages/odepack/dlsode.lisp
+++ b/src/packages/odepack/dlsode.lisp
@@ -68,13 +68,12 @@
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
(defun pend-field (neq time y ydot)
- (setf (fv-ref ydot 0) (fv-ref y 1)
- (fv-ref ydot 1) (- (sin (fv-ref y 0)))))
+ (setf (fv-ref ydot 0) (fv-ref y 1)
+ (fv-ref ydot 1) (- (sin (fv-ref y 0)))))
(defun pend-report (ts y)
(format t "~A ~A ~A ~%" ts (aref y 0) (aref y 1)))
-
(defun pcart-field (neq time y ydot)
(declare (ignore neq time))
(very-quickly
diff --git a/src/level-1/random.lisp b/src/special/random.lisp
similarity index 50%
rename from src/level-1/random.lisp
rename to src/special/random.lisp
index 4bdf231..997d039 100644
--- a/src/level-1/random.lisp
+++ b/src/special/random.lisp
@@ -1,17 +1,18 @@
(in-package #:matlisp)
-(declaim (inline draw-standard-exponential))
-(defun draw-standard-exponential ()
+(declaim (ftype (function () double-float) draw-standard-normal draw-standard-normal-marsaglia draw-standard-exponential))
+
+(definline draw-standard-exponential ()
"Return a random variable from the Exponential(1) distribution, which has density exp(-x)."
+ ;; Adapted from cl-random, originally written by Tamas Papp
;; need 1-random, because there is a small but nonzero chance of getting a 0.
(- (log (- 1d0 (random 1d0)))))
-(declaim (ftype (function () double-float) draw-standard-normal-leva draw-standard-normal-marsaglia))
-(definline draw-standard-normal-leva ()
+(definline draw-standard-normal ()
"Draw a random number from N(0,1)."
;; Method from Leva (1992). This is considered much better/faster than the Box-Muller method.
;; Adapted from cl-random, originally written by Tamas Papp
- ;; This tends to be just as fast as Marsaglia with storage.
+ ;; This seems to be just as fast as Marsaglia with storage.
(very-quickly
(loop
:do (let* ((u (random 1d0))
@@ -43,25 +44,33 @@
(setf prev (* y mult))
(return (* x mult))))))))))
-;;
-(defun randn (dims)
- (let* ((ret (zeros dims 'real-tensor))
- (sto (store ret)))
- (declare (type (simple-array double-float (*)) sto))
- (very-quickly
- (mod-dotimes (idx (dimensions ret))
- :with (linear-sums
- (of-ret (strides ret)))
- :do (setf (aref sto of-ret) (the double-float (draw-standard-normal-leva)))))
- ret))
+;;
+(defmacro fill-tensor (type (func tensor))
+ (using-gensyms (decl (tensor))
+ (with-gensyms (sto ofst)
+ `(let* (,@decl
+ (,sto (store ,tensor)))
+ (declare (type ,type ,tensor)
+ (type ,(store-type type) ,sto))
+ (very-quickly
+ (mod-dotimes (idx (dimensions ,tensor))
+ :with (linear-sums
+ (,ofst (strides ,tensor)))
+ :do (t/store-set ,type ,(etypecase func (symbol `(,func)) (cons func)) ,sto ,ofst)))
+ ,tensor))))
+
+(macrolet ((generate-rand (func clause)
+ (let ((clause (etypecase clause
+ (symbol `(,clause))
+ (cons clause))))
+ `(defun ,func (&optional dims)
+ (if dims
+ (fill-tensor real-tensor (,clause (zeros dims 'real-tensor)))
+ ,clause))))
+ (generate-rands ((&rest args))
+ `(progn
+ ,@(mapcar #'(lambda (x) `(generate-rand ,(car x) ,(cadr x))) args))))
+ (generate-rands ((randn (draw-standard-normal))
+ (rand (random 1d0))
+ (rande (draw-standard-exponential)))))
-(defun rand (dims)
- (let* ((ret (zeros dims 'real-tensor))
- (sto (store ret)))
- (declare (type (simple-array double-float (*)) sto))
- (very-quickly
- (mod-dotimes (idx (dimensions ret))
- :with (linear-sums
- (of-ret (strides ret)))
- :do (setf (aref sto of-ret) (random 1d0))))
- ret))
commit 60f2b8490d5a31c90886c51d081bdb802f5431d6
Author: Akshay Srinivasan <aks...@gm...>
Date: Wed Aug 21 19:50:19 2013 -0700
Made attribute hash-table allocation lazy.
diff --git a/lib-src/gnuplot/gnuplot.lisp b/lib-src/gnuplot/gnuplot.lisp
index 9bfaf56..d34a360 100644
--- a/lib-src/gnuplot/gnuplot.lisp
+++ b/lib-src/gnuplot/gnuplot.lisp
@@ -8,7 +8,7 @@
ccl:run-program
gnuplot-binary nil :input :stream :wait nil :output t))
-(defun plot2d (data &key (color (list "#FF0000")))
+(defun gnuplot-send (str &rest args)
(unless *current-gnuplot-process*
(setf *current-gnuplot-process* (open-gnuplot-stream)))
(let ((stream (#+:sbcl
@@ -16,23 +16,16 @@
#+:ccl
ccl:external-process-input-stream
*current-gnuplot-process*)))
- (with-open-file (s "/tmp/matlisp-gnuplot.out" :direction :output :if-exists :supersede :if-does-not-exist :create)
- (loop :for i :from 0 :below (loop :for x :in data :minimizing (number-of-elements x))
- :do (loop :for x :in data :do (format s "~a " (tensor-ref x i)) :finally (format s "~%"))))
- (format stream "plot '/tmp/matlisp-gnuplot.out' with lines linecolor rgb ~s~%" color)
- (finish-output stream)))
-
-(defun gnuplot-send (str)
- (unless *current-gnuplot-process*
- (setf *current-gnuplot-process* (open-gnuplot-stream)))
- (let ((stream (#+:sbcl
- sb-ext:process-input
- #+:ccl
- ccl:external-process-input-stream
- *current-gnuplot-process*)))
- (format stream "~a~%" str)
+ (apply #'format (append (list stream str) args))
(finish-output stream)))
+(defun plot2d (data &key (lines t) (color (list "#FF0000")))
+ (with-open-file (s "/tmp/matlisp-gnuplot.out" :direction :output :if-exists :supersede :if-does-not-exist :create)
+ (loop :for i :from 0 :below (loop :for x :in data :minimizing (size x))
+ :do (loop :for x :in data :do (format s "~a " (coerce (ref x i) 'single-float)) :finally (format s "~%"))))
+ (if lines
+ (gnuplot-send "plot '/tmp/matlisp-gnuplot.out' with lines linecolor rgb ~s~%" color)
+ (gnuplot-send "plot '/tmp/matlisp-gnuplot.out'~%")))
;; (defclass gnuplot-plot-info ()
;; ((title
diff --git a/src/base/standard-tensor.lisp b/src/base/standard-tensor.lisp
index c34ed6d..426536f 100644
--- a/src/base/standard-tensor.lisp
+++ b/src/base/standard-tensor.lisp
@@ -58,10 +58,22 @@
(store :initarg :store :reader store
:documentation "The actual storage for the tensor.")
;;
- (attributes :initarg :attributes :reader attributes :initform (make-hash-table)
+ (attributes :initarg :attributes :initform nil
:documentation "Place for computable attributes of an object instance."))
(:documentation "Basic tensor class."))
+;;Create hash-table only when necessary
+(definline attributes (x)
+ (declare (type standard-tensor x))
+ (or (slot-value x 'attributes)
+ (let ((htbl (make-hash-table)))
+ (setf (slot-value x 'attributes) htbl)
+ htbl)))
+
+(declaim (ftype (function (standard-tensor) index-store-vector) strides dimensions)
+ (ftype (function (standard-tensor) index-type) head)
+ (ftype (function (standard-tensor) hash-table) attributes))
+
(defmacro memoizing ((tensor name) &rest body)
(declare (type symbol name))
(with-gensyms (tens)
@@ -95,6 +107,7 @@
(declare (type standard-tensor tensor))
(lvec-foldr #'* (the index-store-vector (dimensions tensor))))
+;;
(defgeneric store-size (tensor)
(:documentation "
Syntax
@@ -424,3 +437,14 @@
:strides (very-quickly (vectorify (the index-store-vector nstds) nrank 'index-type))
:store (store tensor)
:parent-tensor tensor)))))))
+
+(definline slice~ (x axis &optional (idx 0))
+ (let ((slst (make-list (rank x) :initial-element '(* * *))))
+ (rplaca (nthcdr axis slst) (list idx '* (1+ idx)))
+ (sub-tensor~ x slst nil)))
+
+(definline row-slice~ (x idx)
+ (slice~ x 0 idx))
+
+(definline col-slice~ (x idx)
+ (slice~ x 1 idx))
diff --git a/src/level-1/map.lisp b/src/level-1/map.lisp
new file mode 100644
index 0000000..3359ea9
--- /dev/null
+++ b/src/level-1/map.lisp
@@ -0,0 +1,72 @@
+(in-package #:matlisp)
+
+(defgeneric mapsor! (func x y)
+ (:documentation "
+ y <- func(x)
+")
+ (:method :before ((func function) (x standard-tensor) (y standard-tensor))
+ (assert (very-quickly (lvec-eq (dimensions x) (dimensions y))) nil 'tensor-dimension-mismatch)))
+
+(defmethod mapsor! ((func function) (x standard-tensor) (y standard-tensor))
+ (let ((clx (class-name (class-of x)))
+ (cly (class-name (class-of y))))
+ (assert (and
+ (member clx *tensor-type-leaves*)
+ (member cly *tensor-type-leaves*))
+ nil 'tensor-abstract-class :tensor-class (list clx cly))
+ (compile-and-eval
+ `(defmethod mapsor! ((func function) (x ,clx) (y ,cly))
+ (let ((sto-x (store x))
+ (sto-y (store y)))
+ (declare (type ,(store-type clx) sto-x)
+ (type ,(store-type cly) sto-y))
+ (very-quickly
+ (mod-dotimes (idx (dimensions x))
+ :with (linear-sums
+ (of-x (strides x))
+ (of-y (strides y)))
+ :do (t/store-set ,cly (funcall func (t/store-ref ,clx sto-x of-x)) sto-y of-y))))
+ y)))
+ (mapsor! func x y))
+
+(defun mapcol (func x &optional (axis 1))
+ (declare (type standard-tensor x))
+ (assert (tensor-matrixp x) nil 'tensor-dimension-mismatch)
+ (let* ((v-x (slice~ x axis))
+ (st-x (aref (the index-store-vector (strides x)) axis))
+ (ret (zeros (aref (the index-store-vector (dimensions x)) axis) (class-of x))))
+ (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
+ :do (progn
+ (setf (ref ret i) (funcall func v-x))
+ (incf (slot-value v-x 'head) st-x)))
+ ret))
+
+(defun tensor-min (x)
+ (let ((min nil))
+ (mod-dotimes (idx (dimensions x))
+ :do (when (or (null min) (< (ref x idx) min))
+ (setf min (ref x idx))))
+ min))
+
+(defun idx-min (v)
+ (let ((min-idx nil)
+ (min nil))
+ (loop :for i :from 0 :below (aref (the index-store-vector (dimensions v)) 0)
+ :do (when (or (null min) (< (ref v i) min))
+ (setf min-idx i
+ min (ref v i))))
+ (values min-idx min)))
+
+(defun mapxis! (func x y &optional (axis 0))
+ (declare (type standard-tensor x y))
+ (multiple-value-bind (view-x view-y) (let ((slst (make-list (rank x) :initial-element '(* * *))))
+ (rplaca (nthcdr axis slst) (list 0 '* 1))
+ (values (sub-tensor~ x slst nil) (sub-tensor~ y slst nil)))
+ (let ((st-x (aref (the index-store-vector (dimensions x)) axis))
+ (st-y (aref (the index-store-vector (dimensions x)) axis)))
+ (loop :for i :from 0 :below (aref (the index-store-vector (dimensions x)) axis)
+ :do (progn
+ ;;May die if you're doing fancy stuff.
+ (funcall func view-x view-y)
+
+
diff --git a/src/level-1/random.lisp b/src/level-1/random.lisp
new file mode 100644
index 0000000..4bdf231
--- /dev/null
+++ b/src/level-1/random.lisp
@@ -0,0 +1,67 @@
+(in-package #:matlisp)
+
+(declaim (inline draw-standard-exponential))
+(defun draw-standard-exponential ()
+ "Return a random variable from the Exponential(1) distribution, which has density exp(-x)."
+ ;; need 1-random, because there is a small but nonzero chance of getting a 0.
+ (- (log (- 1d0 (random 1d0)))))
+
+(declaim (ftype (function () double-float) draw-standard-normal-leva draw-standard-normal-marsaglia))
+(definline draw-standard-normal-leva ()
+ "Draw a random number from N(0,1)."
+ ;; Method from Leva (1992). This is considered much better/faster than the Box-Muller method.
+ ;; Adapted from cl-random, originally written by Tamas Papp
+ ;; This tends to be just as fast as Marsaglia with storage.
+ (very-quickly
+ (loop
+ :do (let* ((u (random 1d0))
+ (v (* 1.7156d0 (- (random 1d0) 0.5d0)))
+ (x (- u 0.449871d0))
+ (y (+ (abs v) 0.386595d0))
+ (q (+ (expt x 2) (* y (- (* 0.19600d0 y) (* 0.25472d0 x))))))
+ (declare (type double-float u v x y q))
+ (unless (and (> q 0.27597d0)
+ (or (> q 0.27846d0)
+ (plusp (+ (expt v 2) (* 4 (expt u 2) (log u))))))
+ (return (/ v u)))))))
+
+;;Not thread safe, obviously
+(let ((prev nil))
+ (defun draw-standard-normal-marsaglia ()
+ (if prev
+ (prog1 prev
+ (setf prev nil))
+ (very-quickly
+ (loop
+ :do (let* ((x (1- (random 2d0)))
+ (y (1- (random 2d0)))
+ (s (+ (* x x) (* y y))))
+ (declare (type double-float x y s))
+ (when (<= s 1d0)
+ (let ((mult (sqrt (/ (* -2 (log s)) s))))
+ (declare (type double-float mult))
+ (setf prev (* y mult))
+ (return (* x mult))))))))))
+
+;;
+(defun randn (dims)
+ (let* ((ret (zeros dims 'real-tensor))
+ (sto (store ret)))
+ (declare (type (simple-array double-float (*)) sto))
+ (very-quickly
+ (mod-dotimes (idx (dimensions ret))
+ :with (linear-sums
+ (of-ret (strides ret)))
+ :do (setf (aref sto of-ret) (the double-float (draw-standard-normal-leva)))))
+ ret))
+
+(defun rand (dims)
+ (let* ((ret (zeros dims 'real-tensor))
+ (sto (store ret)))
+ (declare (type (simple-array double-float (*)) sto))
+ (very-quickly
+ (mod-dotimes (idx (dimensions ret))
+ :with (linear-sums
+ (of-ret (strides ret)))
+ :do (setf (aref sto of-ret) (random 1d0))))
+ ret))
commit 050548a939e93208a4990b6de248e3e39b3caa45
Author: Akshay Srinivasan <aks...@gm...>
Date: Tue Aug 20 01:21:26 2013 -0700
Added sum method. Made changes to dot.
diff --git a/matlisp.asd b/matlisp.asd
index f3c236b..cfb71d3 100644
--- a/matlisp.asd
+++ b/matlisp.asd
@@ -148,7 +148,9 @@
(:file "realimag"
:depends-on ("copy"))
(:file "trans"
- :depends-on ("scal" "copy"))))
+ :depends-on ("scal" "copy"))
+ (:file "sum"
+ :depends-on ("dot" "copy"))))
(:module "matlisp-level-2"
:pathname "level-2"
:depends-on ("matlisp-base" "matlisp-classes" "foreign-core" "matlisp-level-1")
diff --git a/src/level-1/dot.lisp b/src/level-1/dot.lisp
index 354c271..3a81f7e 100644
--- a/src/level-1/dot.lisp
+++ b/src/level-1/dot.lisp
@@ -34,41 +34,48 @@
(if conjp 'zdotc 'zdotu))
;;
-(deft/generic (t/blas-dot #'subtypep) sym (x y &optional conjp))
-(deft/method t/blas-dot (sym blas-numeric-tensor) (x y &optional (conjp t))
- (using-gensyms (decl (x y))
- `(let (,@decl)
- (declare (type ,sym ,x ,y))
- (,(macroexpand-1 `(t/blas-dot-func ,sym ,conjp))
- (aref (the index-store-vector (dimensions ,x)) 0)
- (the ,(store-type sym) (store ,x)) (aref (the index-store-vector (strides ,x)) 0)
- (the ,(store-type sym) (store ,y)) (aref (the index-store-vector (strides ,y)) 0)
- (head ,x) (head ,y)))))
+(deft/generic (t/blas-dot #'subtypep) sym (x y &optional conjp num-y?))
+(deft/method t/blas-dot (sym blas-numeric-tensor) (x y &optional (conjp t) (num-y? nil))
+ (using-gensyms (decl (x y))
+ (with-gensyms (sto)
+ `(let (,@decl
+ ,@(when num-y? `((,sto (t/store-allocator ,sym 1)))))
+ (declare (type ,sym ,x ,@(unless num-y? `(,y)))
+ ,@(when num-y? `((type ,(field-type sym) ,y)
+ (type ,(store-type sym) ,sto))))
+ ,@(when num-y? `((t/store-set ,sym ,y ,sto 0)))
+ (,(macroexpand-1 `(t/blas-dot-func ,sym ,conjp))
+ (aref (the index-store-vector (dimensions ,x)) 0)
+ (the ,(store-type sym) (store ,x)) (aref (the index-store-vector (strides ,x)) 0)
+ (the ,(store-type sym) ,(if num-y? sto `(store ,y))) ,(if num-y? 0 `(aref (the index-store-vector (strides ,y)) 0))
+ (head ,x) ,(if num-y? 0 `(head ,y)))))))
-(deft/generic (t/dot #'subtypep) sym (x y &optional conjp))
-(deft/method t/dot (sym standard-tensor) (x y &optional (conjp t))
+(deft/generic (t/dot #'subtypep) sym (x y &optional conjp num-y?))
+(deft/method t/dot (sym standard-tensor) (x y &optional (conjp t) (num-y? nil))
(using-gensyms (decl (x y))
(with-gensyms (sto-x sto-y of-x of-y stp-x stp-y dot)
`(let (,@decl)
- (declare (type ,sym ,x ,y))
+ (declare (type ,sym ,x ,@(unless num-y? `(,y)))
+ ,@(when num-y? `((type ,(field-type sym) ,y))))
(let ((,sto-x (store ,x))
(,stp-x (aref (the index-store-vector (strides ,x)) 0))
(,of-x (head ,x))
- (,sto-y (store ,y))
- (,stp-y (aref (the index-store-vector (strides ,y)) 0))
- (,of-y (head ,y))
+ ,@(unless num-y?
+ `((,sto-y (store ,y))
+ (,stp-y (aref (the index-store-vector (strides ,y)) 0))
+ (,of-y (head ,y))))
(,dot (t/fid+ ,(field-type sym))))
- (declare (type ,(store-type sym) ,sto-x ,sto-y)
- (type index-type ,stp-x ,stp-y ,of-x ,of-y)
+ (declare (type ,(store-type sym) ,sto-x ,@(unless `(,sto-y)))
+ (type index-type ,stp-x ,of-x ,@(unless num-y? `(,stp-y ,of-y)))
(type ,(field-type sym) ,dot))
(very-quickly
(loop :repeat (aref (the index-store-vector (dimensions ,x)) 0)
:do (setf ,dot (t/f+ ,(field-type sym) ,dot
(t/f* ,(field-type sym)
,(recursive-append (when conjp `(t/fc ,(field-type sym))) `(t/store-ref ,sym ,sto-x ,of-x))
- (t/store-ref ,sym ,sto-y ,of-y)))
+ ,(if num-y? y `(t/store-ref ,sym ,sto-y ,of-y))))
,of-x (+ ,of-x ,stp-x)
- ,of-y (+ ,of-y ,stp-y))))
+ ,@(unless num-y? `(,of-y (+ ,of-y ,stp-y))))))
,dot)))))
;;---------------------------------------------------------------;;
(defgeneric dot (x y &optional conjugate-p)
@@ -133,3 +140,24 @@
(dot x y conjugate-p))
(t
(error "Don't know how to compute the dot product of ~a , ~a." clx cly)))))
+
+(defmethod dot ((x standard-tensor) (y t) &op...
[truncated message content] |
Thanks for helping keep SourceForge clean.
X