Page 1 of 1

Is there any good efficient vector library?

Posted: Fri Mar 27, 2009 6:54 am
by Jasper
Is there any efficient vector library out there? How well is the type (vector &optional element-type size) optimized? (I don't want to sacrifice any performence.) I started making a vector package, but i don't want to duplicate effort. (Well, i had a vector library, but it doesn't satisfy points below.) I can't find any vector libraries, unless they come with other things.(So there is duplication effort in the projects there.

What i expect from a proper vector library:
  • No specifying what type of vector it is anywhere except in the type or in some cases the creation.
  • No storage of integer giving the length when lengths are small. No iterating, expand small loops.
  • v+, v-, v*, v/, inpr, crosspr, lensqr, len. angle-of-v (values if multiple.), v-of-angle (More arguments if multiple.)
  • Preferably (define-vector vector-type maker component-type initform components)-like macro that integrates with the rest. (Same rules for types as with numbers. (* float integer)->float and all. Not shy of integers.
  • Matrix stuff separate package. Might want to call matrix-multiplication the same inpr as the vectors.

Re: Is there any good efficient vector library?

Posted: Fri Mar 27, 2009 7:27 am
by gugamilare
A quick search into http://www.cl-user.net/ gave me the following results:

http://common-lisp.net/project/cl-mathstats/
http://nlisp.info/
http://matlisp.sourceforge.net/
http://common-lisp.net/project/gsll/
http://aleph0.info/spartns/ (for sparse vectors / matrices)

Anyway, it doesn't hurt to program your own library, just for learning and without worrying with performance, but if you want a good one, you can look into those one, and maybe a few more.

Take a look at http://www.cliki.net/ as well.

Re: Is there any good efficient vector library?

Posted: Fri Mar 27, 2009 3:28 pm
by Jasper
Those are a lot of libraries to check. Looking superficially none seem to have what i want.

I made a vector package. It satisfies most of my requirements. It even makes matrices as extensions of macros (inpr being matrix multiplication.) I tried to make the addition/substraction work 'in one go' rather then one by one. I mean doing (+ a b c) as such instead of (+ a (+ b c)) ) but unfortunately that made the method system go crazy, taking ages to load and such, only other way to fix it i see violates my first demand.

Code below is 318 lines total.

vect.lisp Package, basic definitions.

Code: Select all

(defpackage #:vect
  (:use #:common-lisp #:iterate)
  (:export mk-vect define-vect u v w x y z
	   v+ v- v* v/ inpr crosspr lensqr len distsqr dist
	   v-angle angle-v2 angle-v3
	   vect-as-list)
  (:documentation "Vectors/matrices with numbers in them."))

(in-package #:vect)

(defgeneric r* (a b)
  (:documentation "Multiplies two numbers. v* is the macro that users use.\
 There is no r/ it is done via r*."))

(defgeneric r+2 (a b)
  (:documentation "Adds two numbers. v+ is the macro that users use. It 
leads to r+n methods."))
(defgeneric r-2 (a b)
  (:documentation "subtracts two numbers. v+ is the macro that users use.\
 It leads to r+n methods."))

(defgeneric r+3 (a b c) (:documentation "See r+2"))
(defgeneric r+4 (a b c d) (:documentation "See r+2"))
(defgeneric r+5 (a b c d e) (:documentation "See r+2"))
(defgeneric r+6 (a b c d e f) (:documentation "See r+2"))

(defgeneric r-3 (a b c) (:documentation "See r+2"))
(defgeneric r-4 (a b c d) (:documentation "See r+2"))
(defgeneric r-5 (a b c d e) (:documentation "See r+2"))
(defgeneric r-6 (a b c d e f) (:documentation "See r+2"))

(defgeneric mk-vect2 (x y)
  (:documentation "Makes two-dimensional vector, with type inferred from\
 arguments."))

(defgeneric mk-vect3 (x y z)
  (:documentation "Makes three-dimensional vector, with type inferred from\
 arguments."))

(defgeneric mk-vect4 (w x y z)     (:documentation "See mk-vect2"))
(defgeneric mk-vect5 (v w x y z)   (:documentation "See mk-vect2"))
(defgeneric mk-vect6 (u v w x y z) (:documentation "See mk-vect2"))

(defun vect-maker-name (n)
  (case n (2 'mk-vect2) (3 'mk-vect3) (4 'mk-vect4)
	  (5 'mk-vect5) (6 'mk-vect6)))
  
(defmacro mk-vect (&rest args)
  `(,(vect-maker-name (length args)) ,@args))

(defmacro v* (vect &rest factors)
  "Multiplies vect with factors."
  (if (null factors) vect
      `(r* ,vect (v* ,@factors))))

(defmacro v/ (vect &rest dividers)
  "Divides vect with dividers."
  (if (null dividers) vect
      `(r* ,vect (/ (v* ,@dividers)))))

(defmacro v+ (&rest add)
  "Addition macro. Chooses from the r+n functions, chains them if more then\
 6."
  (case (length add)
    (1 (car add))    (2 `(r+2 ,@add)); (3 `(r+3 ,@add))
;    (4 `(r+4 ,@add)) (5 `(r+5 ,@add)) (6 `(r+6 ,@add))
    (t `(v+ (v+2 ,@(subseq add 0 2)) ,@(subseq add 2)))))

(defmacro v- (&rest add)
  "Substraction macro. Chooses from the r-n functions, chains them if\
 more then 6."
  (case (length add)
    (1 (car add))    (2 `(r-2 ,@add)); (3 `(r-3 ,@add))
    (t `(v- (v-2 ,@(subseq add 0 2)) ,@(subseq add 2)))))

;    (4 `(r-4 ,@add)) (5 `(r-5 ,@add)) (6 `(r-6 ,@add))
;    (t `(v- (v-6 ,@(subseq add 0 6)) ,@(subseq add 6)))))


(defgeneric r-inpr (a b) (:documentation "Inproduct between two things.
Used through inpr macro."))

(defgeneric crosspr (a b) (:documentation "Cross product between two\
 things."))

(defmacro inpr (&rest args)
  "Inproduct between things. Usually two."
  `(r-inpr ,(car args) (inpr ,@(cdr args))))

(declaim (inline len lensqr distsqr dist))
(defun lensqr (vect)
  "Length squared of a vector."
  (r-inpr vect vect))

(defun len (vect)
  "Length of a vector."
  (sqrt (lensqr vect)))

(defun distsqr (a b)
  "Distance squared between two vectors."
  (lensqr (v- a b)))
(defun dist (a b)
  "Distance between two vectors."
  (sqrt (distsqr a b)))

(defgeneric v-angle (vect) (:documentation "Angles from a vector."))

(defun angle-v2 (angle)
  "Vector from a angle."
  (mk-vect (cos angle) (sin angle)))

(defun angle-v3 (angle up-angle)
  "Vector from two angles, so three-dimensional.
 (Not the usual illogical polar coordinates, height is sin up-angle.)"
  (mk-vect (* (cos angle) (cos up-angle)) (* (sin angle) (cos up-angle))
	   (sin up-angle)))

(defgeneric vect-as-list (vec) (:documentation "Converts a vector to a list\
 for printing and such."))
vect-number.lisp The vector code treats things internally as 'vectors' so that it can do matrices as well, hence we need to make a few methods to make numbers work.

Code: Select all

(in-package #:vect)

;;Regular number version.
(defmethod r* ((x number) (y number))
  (* x y))

;Addition and substraction.
(eval (let*((vars '(v1 v2 v3 v4 v5 v6))
	     (args (iter (for v in vars)
			 (collect `(,v number)))))
	 `(progn ,@(iter (for n from 2)
			 (for r+ in '(r+2 r+3 r+4 r+5 r+6))
			 (for r- in '(r-2 r-3 r-4 r-5 r-6))
			 (appending `((defmethod ,r+ (,@(subseq args 0 n))
					(+ ,@(subseq vars 0 n)))
				      (defmethod ,r- (,@(subseq args 0 n))
					(+ ,@(subseq vars 0 n)))))))))
;Inproduct.
(defmethod r-inpr ((a number) (b number))
  (* a b))
vect-define.lisp Allows one to make vectors, it connects it to other vectors via mk-vect (which can be overridden)

Code: Select all

(in-package #:vect)

(defun components-from-integer (dimension)
  "Uses standard component names."
  (case dimension (1 '(x)) (2 '(x y)) (3 '(x y z))))

(defmacro do-components (fun maker components &rest args)
  "Macro, does function with the arguments from all the arguments. 
Components is the number of elements, or their names."
  (when (integerp components)
    (setf components (components-from-integer)))
  `(,maker
    ,@(iter (for c in components)
	    (collect
		`(,fun
		  ,@(iter (for a in args)
			 (collect `(,c ,a))))))))

(defmacro make-do-components-method (method fun maker components type
				     arg-cnt)
  "Makes a method from the components."
  (let ((gs (iter (repeat arg-cnt)
		  (collect (gensym)))))
    `(defmethod ,method (,@(iter (for g in gs)
				 (collect `(,g ,type))))
       (do-components ,fun ,maker ,components ,@gs))))

(defun any-eql (eql list)
  "Any of list eql to eql."
  (dolist (el list) (when (eql eql el) (return el))))

(defun all-combinations (list n fun &optional args)
  (dolist (el list)
    (if (= n 1)
      (funcall fun (cons el args))
      (all-combinations list (- n 1) fun (cons el args)))))

(defun r+n-name (n) (case n (2 'r+2) (3 'r+3) (4 'r+4) (5 'r+5) (6 'r+6)))
(defun r-n-name (n) (case n (2 'r-2) (3 'r-3) (4 'r-4) (5 'r-5) (6 'r-6)))

(defvar *vector-types* nil "Types that already exist.")

(defmacro define-vect (name of-type initform components &key maker
			 fun-maker)
  "Defines vector type with components of given type.
Also makes all methods."
  (when (integerp components)
    (setf components (components-from-integer components)))
  
  (unless maker 
    (setf maker (vect-maker-name (length components))))
  
  (when (any-eql name *vector-types*)
    (error "Vector type with this name already exists."))

  (push name *vector-types*)
  
  `(progn
   ;The class.
     (defclass ,name ()
       (,@(iter (for c in components)
		(collect `(,c :initarg ,c :accessor ,c
			      :type ,of-type :initform ,initform))))
       (:documentation ,(format nil "A ~D-dimensional vector of type ~D."
				(length components) of-type)))
   ;Creation.
     (defmethod ,maker (,@(iter (for c in components)
				(collect `(,c ,of-type))))
       (make-instance ',name ,@(iter (for c in components)
				     (appending `(',c ,c)))))
      ,@(when fun-maker
	 `((defun ,fun-maker (,@components)
	     (make-instance ',name ,@(iter (for c in components)
					   (appending `(',c ,c)))))))
   ;Multiplication. (Dividing based on it.)
     (defmethod r* ((vect ,name) (scalar number))
       (mk-vect ,@(iter (for c in components)
			(collect `(* (,c vect) scalar)))))
   ;Adding/substracting, different counts of arguments.
     ,@(iter (for n from 2) (while (<= n 2))
	     (all-combinations *vector-types* n
		 (lambda (tps)
		   (when (any-eql name tps)
		     (let*((vars (iter (repeat n) (collect (gensym))))
			   (args (iter (for v in vars)
				       (for tp in tps)
				       (collect `(,v ,tp)))))
		       (appending
			`((defmethod ,(r+n-name n) (,@args)
			    (do-components v+ mk-vect ,components
					   ,@vars))
			  (defmethod ,(r-n-name n) (,@args)
			    (do-components v- mk-vect ,components
					   ,@vars)))))))))
   ;Other stuff.
     ,@(iter
	(all-combinations *vector-types* 2
	  (lambda (tps)
	    (when (any-eql name tps)
	      (destructuring-bind (tp-a tp-b) tps
		(appending
		 `((defmethod r-inpr ((a ,tp-a) (b ,tp-b)) ;Inproduct
		     (v+ ,@(iter (for c in components)
				 (collect `(r-inpr (,c a) (,c b))))))
		   ,(case (length components)
		      (2 ;Cross product of 2 dimensional is 1 dimensional.
		       (destructuring-bind (x y) components
			 `(defmethod crosspr  ((a ,tp-a) (b ,tp-b))
			    (- (* (,y a) (,x b)) (* (,x a) (,y b))))))
		      (3
		       (destructuring-bind (x y z) components
			 `(defmethod crosspr  ((a ,tp-a) (b ,tp-b))
			    (mk-vect 
			     (- (* (,z a) (,y b)) (* (,y a) (,z b)))
			     (- (* (,x a) (,z b)) (* (,z a) (,x b)))
			     (- (* (,y a) (,x b)) (* (,x a) (,y b))))))))))))))
	(finish))
     ;Listing
     (defmethod vect-as-list ((vec ,name))
       ,(if (any-eql of-type (cdr *vector-types*))
	    `(iter (for el in (do-components identity list ,components vec))
		   (collect (vect-as-list el)))
	    `(do-components identity list ,components vec)))
     ;TODO v-angle
     ))
vect-specific.lisp Some vectors themselves.

Code: Select all

(in-package #:vect)

;;Matrices from vectors.
(define-vect vect2d double-float 0d0 2 :fun-maker mk-vect2d)
(define-vect vect2f single-float 0.0 2 :fun-maker mk-vect2f)
(define-vect vect2i fixnum       0   2 :fun-maker mk-vect2i)

(define-vect vect3d double-float 0d0 3 :fun-maker mk-vect3d)
(define-vect vect3f single-float 0.0 3 :fun-maker mk-vect3f)
(define-vect vect3i fixnum       0   3 :fun-maker mk-vect3i)
Some matrices made from those

Code: Select all

(in-package #:vect)

;;Matrices from vectors.
(define-vect matrix22d vect2d (mk-vect2d 0d0 0d0) 2
	     :fun-maker mk-matrix22d)
(define-vect matrix23d vect3d (mk-vect3d 0d0 0d0 0d0) 2
	       :fun-maker mk-matrix23d)
(define-vect matrix32d vect2d (mk-vect2d 0d0 0d0) 3
	       :fun-maker mk-matrix32d)
(define-vect matrix33d vect3d (mk-vect3d 0d0 0d0 0d0) 3
	       :fun-maker mk-matrix33d)

Re: Is there any good efficient vector library?

Posted: Sat Mar 27, 2010 4:14 pm
by Kohath
Hey, as a note (and I know it's exactly one year from the last post :mrgreen:, so things may have changed since then), I would suggest implementing the r+ and r- functionality with something like this:

Code: Select all

(defun v+ (&rest add)
  (reduce #'r+2 (cdr add) :initial-value (car add)))

(defun v- (&rest add)
  (cond ((cdr add) (reduce #'r-2 (cdr add) :initial-value (car add)))
        (add (negative-of (car add)))
        (t default-value-or-error)))
The function v+ handles the no arguments case by returning nil, but an (or (car add) default-value-or-error) would be one way to fix that if it's a problem.

But if you wanted to do it in one swoop, without calling the 2 argument version, I'd recommend an iterative version, checking the types of each argument and then only if all types are the same (or compatible) would I dispatch to a generic function (which would use the type of one arg to specialize), like this:

Code: Select all

(defmethod same-type-vmult ((prototype (vector integer)) vector-list)
  (loop for vect in vector-list
    do ....))

Re: Is there any good efficient vector library?

Posted: Mon Mar 29, 2010 6:53 am
by Jasper
The newer code can be found in the vect/ directory of the git of this project. Note that i haven't done much at all to that project lately! Infact, i moved some of it to the Gil project so people wouldn't have to download both, and because i can really consider :generic so a lib for other people to use. Some of it is more stuff i like/use and will hopefully come out of a nicely distributed package somewhere later. Consider that lisp-umac thing dead. (Perhaps i should put that on the website)

I would write the code in vect/ rather differently now, although i don't see much way around it. Imo whether it is better to use classes or cl:vectors is pretty much up to the question 'which is optimized better', these days i use vector-vect and hope it will be optimized.

I think having a library(/asd package) focussing solely vanilla vector/matrix operations would be a good idea. vector-vect.lisp in that git repo doesn't do matrices and vectors in one sweep; (for convenience)i'll just copy code that does that below.

Code: Select all

;;Placed in public domain.

(cl:in-package :cl-user)

(defpackage :vector-vect
  (:use :common-lisp)
  (:export vect2d vect3d x y z xyz xyzs with-xyz with-xyzs
	   v+ v- v* v/
	   inpr crosspr lensqr len normalize  distsqr dist
	   angle-v v2-angle rotate-v2 rotate-v3)
  (:documentation "Vectors from sequences.
Hope optimization will do its job well!"))

(in-package :vector-vect)

(deftype vect2d () '(vector number 2))
(deftype vect3d () '(vector number 3))

(defmacro mk-component (n name)
  "Simple utility to make the component functions."
  `(defmacro ,name (vect)
     ,(format nil "~D component of vector. (index ~D)" name n)
     `(aref ,vect ,',n)))

(mk-component 0 x)
(mk-component 1 y)
(mk-component 2 z)

(defmacro with-gensyms ((&rest vars)&body body)
  "Makes you some variables with gensyms output in them."
  `(let (,@(loop for el in vars
	      collect `(,el (gensym))))
     ,@body))

(defmacro xyz (input &body body)
  "Makes settable x,y,z for you, or name-x,y,z"
  (with-gensyms (v)
    (multiple-value-bind (vect name)
	(typecase input
	  (list (values (car input) (cadr input)))
	  (t    (values input nil)))
      (typecase name
	(null ;No name
	 `(let ((,v ,vect))
	    (symbol-macrolet
		((x (aref ,v 0)) (y (aref ,v 1)) (z (aref ,v 2)))
	      ,@body)))
 ;Automatic name. (Note: not sure if providing it is good idea.)
	((eql t)
	 (typecase vect
	   ((symbolp vect) ;Redo it with variable name as symbol.
	    `(xyz (,vect ,vect) ,@body))
	   ((listp vect) ;Use function/macro name to identify.
	    `(xyx (,vect ,(car vect)) ,@body))
	   (t
	    (error "Don't know how to automatically give a name"))))
	((or symbol string) ;Given name.
	 `(symbol-macrolet
	      (,@(loop for c in '(x y z)
		    for k from 0
		    collect `(,(intern(format nil "~D-~D" name c))
			       (aref ,vect ,k))))
	    ,@body))
	(t
	 (error "Didn't know what to do with input."))))))

(defmacro xyzs (inputs &body body)
  "Multiple xyz macros in sequence."
  (if (null (cdr inputs))
    `(xyz ,(car inputs) ,@body)
    `(xyz ,(car inputs) (multi-xyz ,(cdr inputs) ,@body))))

(defmacro with-xyz (input &body body)
  "Synonym for xyz"
  `(xyz ,input ,@body))
(defmacro with-xyzs (inputs &body body)
  "Synonym for xyzs"
  `(xyzs ,inputs ,@body))

(defmacro with-make-array (var dim k &body body)
  `(let ((,var (make-array ,dim :initial-element 0)))
     (dotimes (,k ,dim) ,@body) ;Note abstraction leak on dim.
     ,var))

;;Adding, substracting.
(defun v+ (&rest vects)
  (reduce (lambda (a b)
	    (typecase a
	      (vector
	       (with-make-array out (length a) i
		 (setf (aref out i) (v+ (aref a i) (aref b i)))))
	      (number
	       (+ a b))))
	  vects))

(defun v* (vect &rest factors)
  "Multiplication by _scalars_"
  (let ((factor (reduce #'* factors)))
    (typecase vect
      (vector
       (with-make-array out (length vect) i
	 (setf (aref out i) (* (aref vect i) factor))))
      (number
       (* vect factor)))))

(defun v- (&rest vects)
  (v+ (car vects)
      (v* (apply #'v+ (cdr vects)) -1)))

(defun v/ (vect &rest divide-by)
  "Division by scalars."
  (v* vect (/ (reduce divide-by))))

(declaim (inline v+ v- v* v/))

(defun inpr (&rest vals)
  "Inproducts and matrix multiplication."
  (reduce (lambda (a b)
	    (typecase a
	      (vector
	       (print (list a b))
	       (typecase b
		 (vector
		  (do*((i 0 (+ i 1))
		       (out (inpr (aref a 0) (aref b 0))
			    (v+ out (inpr (aref a i) (aref b i)))))
		      ((>= (+ i 1) (length a)) out)))
		 (number
		  (v* a b))))
	      (number
	       (* a b))))
	  vals))

(defun crosspr (a b)
  "Cross product. Only 2 or 3 dimensional."
  (xyz (a)
    (case (length a)
      (2 (- (* x (y b)) (* y (x b))))
      (3 (vector
	  (- (* y (z b)) (* z (y b)))
	  (- (* z (x b)) (* x (z b)))
	  (- (* x (y b)) (* y (x b))))))))

(defun lensqr (&rest vects)
  "Length squared of vector."
  (loop for v in vects sum (inpr v v)))
(defun len (&rest vects)
  "Length of vector."
  (sqrt (apply #'lensqr v)))

(defun normalize (v)
  (v/ v (len v)))

(declaim (inline inpr crosspr lensqr len normalize))

;;A little bit arbitrary to have this one generic.. the others can be seen abstractly too..
(defgeneric distsqr (a b)
  (:documentation "Distance squared between two things."))
(defmethod distsqr ((a vector) (b vector))
  (lensqr (v- a b)))
(defgeneric dist (a b)
  (:documentation "Distance between two things."))
(defmethod dist (a b)
  (sqrt (distsqr a b)))

(defun angle-v-raw (angles &key pre)
  (if (null angles)
    `((* ,@pre))
    `((* ,@pre (cos ,(car angles)))
      ,@(angle-v-raw (cdr angles) :pre `(,@pre (sin ,(car angles)))))))

(defmacro angle-v (&rest angles)
  "Produces an n+1 dimensional vector from n angles."
  (let ((angle-vars (loop for el in angles
			 collect (gensym))))
    `(let (,@(loop for el in angles
		   for gs in angle-vars
		collect `(,gs ,el)))
       `(vector ,@(angle-v-raw angle-vars)))))

(defun v2-angle (v)
  "Gets angle from 2d vector."
  (let ((len (len v)))
    (xyz v
      (if (= len 0) 0
	(* (acos (/ x len)) (if (< y 0) -1 1))))))

(defun rotate-v2 (v a &key (cos-a (cos a)) (sin-a (sin a)))
  "Rotates a 2d vector."
  (xyz v
    (vector (+ (* y sin-a) (* x cos-a))
	    (- (* y cos-a) (* x sin-a)))))

(defun rotate-v3 (v n a &key (cos-a (cos a)) (sin-a (sin a)))
  "Rotate 3d vector round something"
  (declare (type vect3d v n)
	   (type number a cos-a sin-a))
  (let*((e1 (xyz (n)
	      (normalize
	       (cond
		 ((and (= x 0) (= y 0))
		  (return-from rotate-v3
		    (xyz (v)
		      (vector (+ (* y sin-a) (* x cos-a))
			      (- (* y cos-a) (* x sin-a))
			      z)))) 
		 ((= x 0)
		  (vector 0d0 (/ -1 y) (/ z)))
		 ((= y 0)
		  (vector (/ -1 x) 0d0 (/ z)))
		 (t
		  (vector (/ x) (/ -1 y) 0d0))))))
	(e2 (normalize (crosspr e1 n)))
	(x (inpr e1 v))
	(y (inpr e2 v)))
    (declare (type vect3d e1 e2)
	     (type number x y))
    (v+ (v* e1 (+ (* y sin-a) (* x cos-a)))
	(v* e2 (- (* y cos-a) (* x sin-a)))
	(v* n  (inpr (normalize n) v)))))

(defun rotate-matrix-2 (a &key (cos-a (cos a)) (sin-a (sin a)))
  "Make the matrix to rotate"
  (declare (type number a cos-a cos-b))
  (vector (vector cos-a     sin-a)
	  (vector (- sin-a) cos-a)))

;TODO defun rotate-matrix-3, matrix stuff etcetera.