3 ;;; Copyright (c) 2004-2006 Cyrus Harmon (ch-lisp@bobobeach.com)
4 ;;; All rights reserved.
6 ;;; Redistribution and use in source and binary forms, with or without
7 ;;; modification, are permitted provided that the following conditions
10 ;;; * Redistributions of source code must retain the above copyright
11 ;;; notice, this list of conditions and the following disclaimer.
13 ;;; * Redistributions in binary form must reproduce the above
14 ;;; copyright notice, this list of conditions and the following
15 ;;; disclaimer in the documentation and/or other materials
16 ;;; provided with the distribution.
18 ;;; THIS SOFTWARE IS PROVIDED BY THE AUTHOR 'AS IS' AND ANY EXPRESSED
19 ;;; OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
20 ;;; WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 ;;; ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
22 ;;; DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 ;;; DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
24 ;;; GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25 ;;; INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
26 ;;; WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
27 ;;; NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
28 ;;; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 (:documentation
"Returns a list containg the number of
38 elments in each dimension of the matrix."))
41 (:documentation
"Returns the number of rows in the matrix."))
44 (:documentation
"Returns the number of columns in the matrix."))
46 (defgeneric val
(m i j
)
47 (:documentation
"Returns the value of the element in the ith row of
48 the jth column of the matrix m."))
50 ;;; This is totally bogys. Why do we have both mref and val? I'm
51 ;;; assuming that val should go away!
52 (defgeneric mref
(m &rest indices
)
53 (:documentation
"Returns the value of the element in the ith row of
54 the jth column of the matrix m."))
56 (defgeneric (setf mref
) (v m
&rest indices
)
57 (:documentation
"Set the value of the specified element at row row
58 and col col of matrix m to be v."))
60 (defgeneric row-major-mref
(m index
))
62 (defgeneric (setf row-major-mref
) (v m index
))
64 (defgeneric move-element
(m i1 j1 n i2 j2
)
65 (:documentation
"Copy the contents of the element at row i1, column
66 j1, in matrix m to the element at row i2, column j2, in matrix n."))
68 (defgeneric set-val
(m i j v
&key coerce
)
69 (:documentation
"Sets the value of the element at row i, column j of
72 (defgeneric set-val
(m i j v
&key coerce
))
74 ;;;; level-1 interface
76 ;;; arithmetic functions
78 (defgeneric m
+ (&rest matrices
)
79 (:documentation
"Element-wise addition of matrices. All matrices
80 must be of the same dimensions. Returns a matrix of the same size as
81 each matrix in matrices, whoses elements are the some of the
82 corresponding elements from each matrix in matrices."))
84 (defgeneric m-
(&rest matrices
)
85 (:documentation
"When passed a single matrix, returns a new
86 matrix of the same dimensions as matrix whose elemnts are the
87 result of perforing a unary minus (or sign inversion) on each
88 element in the matrix. When passed more than one matrix,
89 performs element-wise subtraction of each matrix after the
90 first from the first matrix and returns a new matrix containing
93 (defgeneric m
* (&rest matrices
)
94 (:documentation
"General purpose matrix multiplication
95 operator. If one argument is supplied, returns that matrix. If
96 two arguments are supplied, if both are matrices, performs a
97 matrix multiplication equivalent to multiplying the first
98 matrix by the second matrix. if one argument is a matrix and
99 the other a number, m* scales the matrix by the numeric
100 argument. If more than two arguments are supplied, the first
101 two arguments are treated as in the two argument case and the
102 results are then multiplied by the remaining arguments such
103 that (m* m1 m2 m3) == (m* (m* m1 m2) m3)."))
105 (defgeneric m.
* (&rest matrices
)
106 (:documentation
"Hadamard multiplication operator. Performs an
107 element-wise multiplication of the elements in each matrix."))
109 ;;; logical functions
111 (defgeneric mlognot
(m &key in-place
)
112 (:documentation
"Performs element-wise logical negation of the
113 matrix m. If in-place is nil, returns a new matrix with the
114 resulting values, otherwise, destructively modifies matrix
117 (defgeneric mlognot-range
(m startr endr startc endc
&key in-place
))
119 (defgeneric mlogand-range
(m1 m2 startr endr startc endc
&key in-place
))
120 (defgeneric mlogand
(m1 m2
&key in-place
))
122 (defgeneric mlogior-range
(m1 m2 startr endr startc endc
&key in-place
))
123 (defgeneric mlogior
(m1 m2
&key in-place
))
125 (defgeneric mlogxor-range
(m1 m2 startr endr startc endc
&key in-place
))
126 (defgeneric mlogxor
(m1 m2
&key in-place
))
128 ;;; FIXME!! fix bitnor
129 (defgeneric mbitnor-range
(m1 m2 startr endr startc endc
))
130 (defgeneric mbitnor
(m1 m2
))
132 ;;; extrema (min and max)
134 (defgeneric min-range
(m startr endr startc endc
))
135 (defgeneric max-range
(m startr endr startc endc
))
136 (defgeneric sum-range
(m startr endr startc endc
))
137 (defgeneric min-val
(m))
138 (defgeneric max-val
(m))
140 ;;; exponential functions
142 (defgeneric mat-square
(u))
143 (defgeneric mat-square
! (u))
144 (defgeneric mat-sqrt
(u))
145 (defgeneric mat-sqrt
! (u))
150 (defgeneric sum-cols
(m &key matrix-class
))
151 (defgeneric sum-rows
(m &key matrix-class
))
152 (defgeneric sum-square-range
(m startr endr startc endc
))
153 (defgeneric sum-square
(m))
157 (defgeneric mean-range
(m startr endr startc endc
))
158 (defgeneric mean
(m))
159 (defgeneric variance-range
(m startr endr startc endc
))
160 (defgeneric variance
(m))
161 (defgeneric sample-variance-range
(m startr endr startc endc
))
162 (defgeneric sample-variance
(m))
167 (defgeneric normalize
(u &key normin normax copy
))
168 (defgeneric norm-0-255
(u &key copy
))
169 (defgeneric norm-0-1
(u &key copy
))
171 ;;;; level-0 interface
173 (defgeneric matrix-move-range-2d
(m n startr1 endr1 startc1 endc1
174 startr2 endr2 startc2 endc2
))
176 (defgeneric matrix-move-range-2d-constrain
177 (m n startr1 endr1 startc1 endc1
178 startr2 endr2 startc2 endc2
))
180 (defgeneric matrix-move
(m n
&key constrain
))
182 (defgeneric mat-mult-3-block
(m n p
))
185 ;;; log generic functions
187 (defgeneric mlog-range
(m startr endr startc endc
&optional base
))
189 (defgeneric mlog-range
! (m startr endr startc endc
&optional base
))
191 (defgeneric mlog
(m &optional base
))
193 (defgeneric mlog
! (m &optional base
))
197 (defgeneric mat-hprod-range
(m n startr endr startc endc
))
199 (defgeneric mat-hprod
(m n
))
201 (defgeneric mat-hprod-range
! (m n startr endr startc endc
))
203 (defgeneric mat-hprod
! (m n
))
207 (defgeneric mat-add
(a b
&key in-place
))
208 (defgeneric mat-add-range
(m n startr endr startc endc
&key in-place
))
212 (defgeneric mabs
(u))
214 (defgeneric mabs-range
(m startr endr startc endc
))
216 (defgeneric mabs-range
! (m startr endr startc endc
))
220 (defgeneric %transfrom-matrix
(m n xfrm
&key background interpolation
))
222 ;;; discrete convolution
224 (defgeneric %discrete-convolve
(u v z
&key norm-v
))
226 (defgeneric %separable-discrete-convolve
(m h1 h2 z1 z2
&key truncate norm-v matrix-class
))
230 ;;; Miscellaneous class utilities
233 (defun compute-class-precedence-list (class)
234 (ccl:class-precedence-list class
))
236 (defun subclassp (c1 c2
)
237 (subtypep (class-name c1
) (class-name c2
)))
239 (defgeneric matrix-precedence-list
(c)
241 (let ((mpl (compute-class-precedence-list c
))
242 (mc (find-class 'matrix
)))
243 (mapcan #'(lambda (x)
244 (if (subclassp x mc
) (when x
(list x
))))
248 (defgeneric closest-common-matrix-class
(m1 &rest mr
)
249 (:method
((m1 matrix
) &rest mr
)
250 (car (apply #'ch-util
:closest-common-ancestor
251 (mapcar #'(lambda (x) (matrix-precedence-list (class-of x
)))
254 ;;; gory implementation details follow
256 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
260 (defgeneric fit
(m val
))
261 (defmethod fit ((m matrix
) val
)
265 (defgeneric fit-value
(val m
))
266 (defmethod fit-value (val (m matrix
))
270 (defmethod dim ((m matrix
)) (matrix-dimensions m
))
272 (defmethod rows ((m matrix
))
273 (let ((vals (matrix-vals m
)))
275 (> (length (array-dimensions vals
)) 0))
276 (array-dimension vals
0)
279 (defmethod cols ((m matrix
))
280 (let ((vals (matrix-vals m
)))
282 (> (length (array-dimensions vals
)) 1))
283 (array-dimension vals
1)
286 (declaim (inline set-val
))
287 (defmethod val ((m matrix
) i j
) (aref (matrix-vals m
) i j
))
290 ;;; rvref treats the matrix as a row-vector (that is a
291 ;;; 1 x n matrix). we should throw an error if this is not the case.
292 ;;; TODO: add rvref methods for row-vector, column-vector, and scalar
293 ;;; classes, erroring as appropriate.
294 (defgeneric rvref
(rv i
))
295 (defmethod rvref ((rv matrix
) i
) (aref (matrix-vals rv
) 1 i
))
297 ;;; cvref treats the matrix as a row-vector (that is an
298 ;;; n x 1 matrix). we should throw an error if this is not the case.
299 ;;; TODO: add cvref methods for row-vector, column-vector, and scalar
300 ;;; classes, erroring as appropriate.
301 (defgeneric cvref
(cv i
))
302 (defmethod cvref ((cv matrix
) i
) (aref (matrix-vals cv
) i
1))
305 (defmethod move-element ((m matrix
) i1 j1
(n matrix
) i2 j2
)
306 (setf (mref n i2 j2
) (mref m i1 j1
)))
308 (defmethod set-val ((m matrix
) i j v
&key
(coerce t
))
309 (setf (aref (matrix-vals m
) i j
)
311 (coerce v
(element-type (class-of m
)))
314 (defmethod set-val-fit ((m matrix
) i j v
&key
(truncate nil
))
315 (setf (mref m i j
) (if truncate
(truncate v
) v
)))
318 (defmethod transpose ((m matrix
))
319 (destructuring-bind (rows cols
) (dim m
)
320 (declare (type fixnum rows cols
))
321 (let ((tr (make-instance (class-of m
) :rows cols
:cols rows
)))
323 (declare (type fixnum i
))
325 (declare (type fixnum j
))
326 (set-val tr j i
(val m i j
))))
329 (defmethod transpose ((m double-float-matrix
))
330 (destructuring-bind (rows cols
) (dim m
)
331 (declare (type fixnum rows cols
))
332 (let ((tr (make-instance 'double-float-matrix
:rows cols
:cols rows
)))
333 (with-typed-mref (m double-float
)
334 (with-typed-mref (tr double-float
)
337 (setf (mref tr j i
) (mref m i j
))))))
340 (defgeneric mat-mult
(a b
))
341 (defmethod mat-mult ((a matrix
) (b matrix
))
342 (destructuring-bind (ar ac
) (dim a
)
343 (destructuring-bind (br bc
) (dim b
)
345 (let* ((c (make-instance (class-of a
) :rows ar
:cols bc
)))
350 (incf v
(* (val a i r
) (val b r j
))))
353 (error 'matrix-argument-error
355 "Incompatible matrix dimensions in mat-mult (~S x ~S) * (~S x ~S)."
356 :format-arguments
(list ar ac br bc
))))))
358 (declaim (inline inferior
))
359 (defun inferior (a b
)
362 (declaim (inline superior
))
363 (defun superior (a b
)
366 (declaim (inline constrain
))
367 (defun constrain (min val max
)
369 (inferior (superior min val
) max
)
372 ;;; FIXME - figure out what to do about truncate
373 (defgeneric mat-copy-into
(a c
&key truncate constrain
))
374 (defmethod mat-copy-into ((a matrix
) (c matrix
) &key truncate constrain
)
375 (declare (ignore truncate constrain
))
376 (destructuring-bind (m n
) (dim a
)
378 (declare (type fixnum i
))
380 (declare (type fixnum j
))
381 (move-element a i j c i j
)))
384 (defgeneric mat-copy-proto-dim
(a m n
&key initial-element
))
385 (defmethod mat-copy-proto-dim ((a matrix
) (m fixnum
) (n fixnum
)
386 &key
(initial-element nil initial-element-supplied-p
))
387 (apply #'make-instance
(class-of a
) :rows m
:cols n
388 (when initial-element-supplied-p
389 (list :initial-element initial-element
))))
392 (defgeneric mat-copy-proto
(a))
393 (defmethod mat-copy-proto ((a matrix
))
394 (destructuring-bind (m n
) (dim a
)
395 (make-instance (class-of a
) :rows m
:cols n
)))
397 (defgeneric mat-copy
(a &rest args
))
398 (defmethod mat-copy ((a matrix
) &rest args
)
399 (let ((c (mat-copy-proto a
)))
400 (apply #'mat-copy-into a c args
)
403 (defgeneric mat-scalar-op
(a b op
))
404 (defmethod mat-scalar-op ((a matrix
) (b matrix
) op
)
405 (and (equal (dim a
) (dim b
))
406 (destructuring-bind (m n
) (dim a
)
407 (let ((c (mat-copy a
)))
410 (set-val c i j
(funcall op
(val a i j
) (val b i j
)))))))))
412 (defmethod mat-scalar-op ((a number
) (b matrix
) op
)
413 (destructuring-bind (m n
) (dim b
)
414 (let ((c (mat-copy b
)))
417 (set-val c i j
(funcall op a
(val b i j
))))))))
419 (defmethod mat-scalar-op ((a matrix
) (b number
) op
)
420 (destructuring-bind (m n
) (dim a
)
421 (let ((c (mat-copy a
)))
424 (set-val c i j
(funcall op
(val a i j
) b
)))))))
426 (defgeneric mat-subtr
(a b
&key in-place result-type
))
428 (defgeneric swap-rows
(a k l
))
429 (defmethod swap-rows ((a matrix
) k l
)
433 (let ((h (val a k j
)))
434 (set-val a k j
(val a l j
))
435 (set-val a l j h
)))))
437 (defgeneric swap-cols
(a k l
))
438 (defmethod swap-cols ((a matrix
) k l
)
442 (let ((h (val a i k
)))
443 (set-val a i k
(val a i l
))
444 (set-val a i l h
)))))
446 (defgeneric map-col
(a k f
))
447 (defmethod map-col ((a matrix
) k f
)
448 (destructuring-bind (m n
) (dim a
)
451 (set-val a i k
(funcall f
(val a i k
))))))
453 (defgeneric map-row
(a k f
))
454 (defmethod map-row ((a matrix
) k f
)
455 (destructuring-bind (m n
) (dim a
)
458 (set-val a k j
(funcall f
(val a k j
))))))
460 (macrolet ((frob (name op
)
462 ;;; define the row op
463 (defgeneric ,(ch-util:make-intern
(concatenate 'string
"scalar-" name
"-row"))
465 (defmethod ,(ch-util:make-intern
(concatenate 'string
"scalar-" name
"-row"))
467 (map-row a k
#'(lambda (x) (apply ,op
(list x q
))))
469 ;;; define the column op
470 (defgeneric ,(ch-util:make-intern
(concatenate 'string
"scalar-" name
"-col"))
472 (defmethod ,(ch-util:make-intern
(concatenate 'string
"scalar-" name
"-col"))
474 (map-col a k
#'(lambda (x) (apply ,op
(list x q
))))
478 (frob "double-float-divide" #'/)
479 (frob "single-float-divide" #'/))
481 (defgeneric scalar-mult
(m q
))
482 (defmethod scalar-mult ((m matrix
) q
)
483 (dotimes (i (first (dim m
)) m
)
484 (scalar-mult-row m i q
))
487 (defgeneric scalar-mult-copy
(m q
))
488 (defmethod scalar-mult-copy ((a matrix
) q
)
489 (let ((m (mat-copy a
)))
492 (defgeneric scalar-divide
(a q
))
493 (defmethod scalar-divide ((a matrix
) q
)
494 (dotimes (i (first (dim a
)) a
)
495 (scalar-divide-row a i q
))
498 (defgeneric scalar-divide-copy
(a q
))
499 (defmethod scalar-divide-copy ((a matrix
) q
)
500 (let ((m (mat-copy a
)))
501 (scalar-divide m q
)))
503 (defgeneric zero-matrix
(j k
&key matrix-class
))
504 (defmethod zero-matrix ((j fixnum
) (k fixnum
) &key
(matrix-class 'matrix
))
505 (make-instance matrix-class
:rows j
:cols k
))
507 (defgeneric identity-matrix
(k &key matrix-class
))
508 (defmethod identity-matrix ((k fixnum
) &key
(matrix-class 'matrix
))
509 (let* ((a (zero-matrix k k
:matrix-class matrix-class
))
510 (one (coerce 1 (element-type (class-of a
)))))
512 (set-val a i i one
))))
514 (defgeneric concat-matrix-cols
(a b
&key matrix-type
))
515 (defmethod concat-matrix-cols ((a matrix
) (b matrix
) &key matrix-type
)
516 (let ((da (dim a
)) (db (dim b
)))
518 ((equal (first da
) (first db
))
519 (let ((c (make-instance (if matrix-type
521 (class-of a
)) :rows
(first da
) :cols
(+ (second da
) (second db
)))))
527 (set-val c i j
(val a i j
) :coerce t
)))
530 (set-val c i
(+ j n
) (val b i j
) :coerce t
))))
534 (defgeneric subset-matrix-cols
(a x y
&key
(matrix-type)))
535 (defmethod subset-matrix-cols ((a matrix
) x y
&key
(matrix-type))
539 (let* ((m (first da
))
541 (c (make-instance (if matrix-type
543 (class-of a
)) :rows
(first da
) :cols
(- y x
))))
546 (set-val c i j
(val a i
(+ j x
)))))
550 ;;; this should be optimized via a type-specific method
551 ;;; for each matrix class
552 (defgeneric get-first-non-zero-row-in-col
(a j
&optional start
))
553 (defmethod get-first-non-zero-row-in-col ((a matrix
) j
&optional
(start 0))
554 (let ((n (first (dim a
))))
555 (do ((i start
(+ i
1)))
557 (not (= (val a i j
) 0)))
562 ;;; this should be optimized via a type-specific method
563 ;;; for each matrix class
564 (defgeneric invert-matrix
(a))
565 (defmethod invert-matrix ((a matrix
))
566 (let* ((n (second (dim a
)))
567 (c (concat-matrix-cols a
(identity-matrix n
) :matrix-type
'double-float-matrix
)))
570 ((or (= y n
) (not c
)))
571 (let ((z (get-first-non-zero-row-in-col c y y
)))
573 ((not z
) (setf c nil
))
577 (scalar-divide-row c y
(val c y y
))
582 (let ((k (val c i y
)))
584 (set-val c i j
(+ (val c i j
)
585 (* (- k
) (val c y j
))))))))))))
587 (subset-matrix-cols c n
(+ n n
) :matrix-type
'double-float-matrix
))))
589 (defgeneric transpose-matrix
(a))
590 (defmethod transpose-matrix ((a matrix
))
594 (c (make-instance (class-of a
) :rows n
:cols m
)))
597 (set-val c j i
(val a i j
))))
600 (defgeneric add-row
(m &key values initial-element
)
601 (:method
((m matrix
) &key values initial-element
)
604 (if (null initial-element
)
605 (setf initial-element
(initial-element m
)))
607 (setf (matrix-vals m
)
608 (adjust-array (matrix-vals m
) (list (+ 1 (first d
)) (second d
))
609 :initial-element initial-element
))
615 (set-val m
(first d
) i
(car l
))))))
616 (error 'matrix-error
:message
"Tried to add-row to non-adjustable matrix ~A" m
))))
618 (defgeneric add-col
(m &key values initial-element
)
619 (:method
((m matrix
) &key values initial-element
)
622 (if (null initial-element
)
623 (setf initial-element
(initial-element m
)))
625 (setf (matrix-vals m
)
626 (adjust-array (matrix-vals m
) (list (first d
) (+ 1 (second d
)))
627 :initial-element initial-element
))
633 (set-val m i
(second d
) (car l
))))))
634 (error 'matrix-error
:message
"Tried to add-col to non-adjustable matrix ~A" m
))))
636 (defgeneric reshape
(m rows cols
&key initial-element
)
637 (:method
((m matrix
) rows cols
&key initial-element
)
640 (if (null initial-element
)
641 (setf initial-element
(initial-element m
)))
642 (setf (matrix-vals m
)
643 (adjust-array (matrix-vals m
) (list rows cols
)
644 :initial-element initial-element
)))
645 (error 'matrix-error
:message
"Tried to reshape non-adjustable matrix ~A" m
))))
647 (defgeneric horzcat
(m1 &rest mr
)
648 (:method
((m1 matrix
) &rest mr
)
649 (let ((rows (apply #'max
(mapcar #'rows
(cons m1 mr
))))
650 (cols (apply #'+ (mapcar #'cols
(cons m1 mr
)))))
651 (let ((z (make-instance (apply #'closest-common-matrix-class m1 mr
)
655 (dolist (x (cons m1 mr
))
656 (map-set-range z
0 (- (rows x
) 1) coff
(+ coff
(- (cols x
) 1))
657 #'(lambda (v i j
) (declare (ignore v
)) (val x i
(- j coff
))))
662 (defgeneric vertcat
(m1 &rest mr
)
663 (:method
((m1 matrix
) &rest mr
)
664 (let ((rows (apply #'+ (mapcar #'rows
(cons m1 mr
))))
665 (cols (apply #'max
(mapcar #'cols
(cons m1 mr
)))))
666 (let ((z (make-instance (apply #'closest-common-matrix-class m1 mr
)
670 (dolist (x (cons m1 mr
))
671 (map-set-range z roff
(+ roff
(- (rows x
) 1)) 0 (- (cols x
) 1)
672 #'(lambda (v i j
) (declare (ignore v
)) (val x
(- i roff
) j
)))
677 (defgeneric pad-matrix
(m))
678 (defmethod pad-matrix ((m matrix
))
680 ((> (rows m
) (cols m
))
681 (let ((delta (- (rows m
) (cols m
))))
682 (horzcat (make-instance (class-of m
)
684 :cols
(ceiling (/ delta
2)))
686 (make-instance (class-of m
)
688 :cols
(floor (/ delta
2))))))
689 ((> (cols m
) (rows m
))
690 (let ((delta (- (cols m
) (rows m
))))
691 (vertcat (make-instance (class-of m
)
692 :rows
(ceiling (/ delta
2))
695 (make-instance (class-of m
)
696 :rows
(floor (/ delta
2))
699 (defgeneric set-row
(m r values
))
700 (defmethod set-row ((m matrix
) r values
)
705 (set-val m r i
(car l
))))
707 (defgeneric set-col
(m c values
))
709 (defmethod set-col ((m matrix
) c values
)
714 (set-val m i c
(car l
))))
716 (defgeneric get-row-list
(m r
&optional start
))
718 (defmethod get-row-list ((m matrix
) r
&optional
(start 0))
720 ((< start
(second (dim m
)))
721 (cons (val m r start
)
722 (get-row-list m r
(+ 1 start
))))
725 (defgeneric get-col-list
(m c
&optional start
))
727 (defmethod get-col-list ((m matrix
) c
&optional
(start 0))
729 ((< start
(first (dim m
)))
730 (cons (val m start c
)
731 (get-col-list m c
(+ 1 start
))))
734 (defgeneric map-set-val
(a f
))
736 (defmethod map-set-val ((m matrix
) f
)
737 (loop for i from
0 below
(matrix-total-size m
)
738 do
(setf (row-major-mref m i
)
739 (funcall f
(row-major-mref m i
))))
742 ;;; FIXME THIS IS BROKEN!
743 (defgeneric map-set-val-fit
(a f
&key truncate
))
744 (defmethod map-set-val-fit ((a matrix
) f
&key
(truncate t
))
745 (destructuring-bind (m n
) (dim a
)
746 (declare (dynamic-extent m n
) (fixnum m n
))
748 (declare (dynamic-extent i
) (fixnum i
))
750 (declare (dynamic-extent j
) (fixnum j
))
751 (set-val-fit a i j
(funcall f
(val a i j
)) :truncate truncate
))))
754 (defgeneric map-set-val-copy
(a f
))
755 (defmethod map-set-val-copy ((a matrix
) f
)
756 (destructuring-bind (ar ac
) (dim a
)
757 (declare (dynamic-extent ar ac
) (fixnum ar ac
))
758 (let* ((b (mat-copy-proto a
)))
760 (declare (dynamic-extent i
) (fixnum i
))
762 (declare (dynamic-extent j
) (fixnum j
))
763 (set-val b i j
(funcall f
(val a i j
)))))
766 (defgeneric map-range
(a startr endr startc endc f
))
767 (defmethod map-range ((a matrix
)
773 (declare (dynamic-extent startr endr startc endc
)
774 (fixnum startr endr startc endc
))
775 (do ((i startr
(1+ i
)))
777 (declare (dynamic-extent i
) (type fixnum i
))
778 (do ((j startc
(1+ j
)))
780 (declare (dynamic-extent j
) (type fixnum j
))
781 (funcall f
(val a i j
) i j
))))
783 (defgeneric map-set-range
(a startr endr startc endc func
))
784 (defmethod map-set-range ((a matrix
)
790 (declare (dynamic-extent startr endr startc endc
)
791 (fixnum startr endr startc endc
))
792 (do ((i startr
(1+ i
)))
794 (declare (dynamic-extent i
) (fixnum i
))
795 (do ((j startc
(1+ j
)))
797 (declare (dynamic-extent j
) (fixnum j
))
798 (set-val a i j
(funcall f
(val a i j
) i j
)))))
800 (defgeneric random-matrix
(rows cols
&key matrix-class limit
)
801 (:documentation
"Create a matrix of type <matrix-class> having
802 <rows> rows and <cols> cols. The values of the matrix will be
803 random numbers of the appropriate type between 0 and <limit>."))
805 (defmethod random-matrix ((rows fixnum
) (cols fixnum
) &key
806 (matrix-class 'matrix
)
808 (declare (dynamic-extent rows cols
)
810 (let ((a (make-instance matrix-class
:rows rows
:cols cols
)))
811 (map-set-val a
#'(lambda (x) (declare (ignore x
)) (random limit
)))
815 (defun count-range (startr endr startc endc
)
816 (* (1+ (- endr startr
)) (1+ (- endc startc
))))
818 (defun double-float-divide (&rest args
)
819 (apply #'/ (mapcar #'(lambda (x) (coerce x
'double-float
)) args
)))
821 (defgeneric subset-matrix
(u startr endr startc endc
))
822 (defmethod subset-matrix ((u matrix
) startr endr startc endc
)
823 (destructuring-bind (ur uc
) (dim u
)
825 ((and (<= startr endr ur
) (<= startc endc uc
))
826 (let* ((m (1+ (- endr startr
)))
827 (n (1+ (- endc startc
)))
828 (c (mat-copy-proto-dim u m n
)))
831 (set-val c i j
(val u
(+ i startr
) (+ j startc
)))))
835 (defgeneric map-matrix-copy
(a f
&key matrix-class
))
836 (defmethod map-matrix-copy ((a matrix
) f
&key
(matrix-class (class-of a
)))
837 (destructuring-bind (m n
) (dim a
)
838 (let* ((b (make-instance matrix-class
:rows m
:cols n
)))
841 (set-val b i j
(funcall f a i j
))))
844 (defgeneric array-
>matrix
(a &key matrix-class
))
845 (defmethod array->matrix
((a array
) &key
(matrix-class 'matrix
))
846 (let ((d (array-dimensions a
))
847 (element-type (element-type (find-class matrix-class
)))
849 (cond ((= (length d
) 2)
850 (destructuring-bind (ar ac
) d
852 #|
((and (= ar
1) (= ac
1))
853 (setf m
(scalar (aref a
0 0))))
855 (setf m
(array->col-vector a
)))
857 (setf m
(array->row-vector a
)))
860 (setf m
(make-instance matrix-class
:rows ar
:cols ac
))
864 (set-val m i j
(aref a i j
) :coerce t
)
867 (setf m
(make-instance matrix-class
:dimensions d
))
868 (loop for i from
0 below
(matrix-total-size m
)
869 do
(setf (row-major-mref m i
)
875 (defgeneric matrix-
>list
(m))
876 (defmethod matrix->list
((m matrix
))
877 (destructuring-bind (mr mc
) (dim m
)
879 append
(loop for j below mc
880 collect
(mref m i j
)))))
882 (defgeneric mat-trim
(m k
))
883 (defmethod mat-trim ((m matrix
) k
)
884 (destructuring-bind (mr mc
) (dim m
)
885 (subset-matrix m k
(- mr k
1) k
(- mc k
1))))