Thanks to all who helped with my attempt at producing array numerical operators for lisp. I now
have some code that runs! This code implements 'array syntax' like those in numerical languages
like yorick, matlab, FORTRAN 95. The syntax is to put a decimal point in front of the operator,
so '.+' is array +, '.sin' is array sin etc. The operators .+ and .* work with any mixture and
number of array and scaler arguments. The operators .- and ./ work with only 2 arguments.
Although it works, it's horribly inefficient. I'll have a go at making it run faster next, and
other improvements as I learn lisp. (This is my first lisp program...) I'll start by choosing
the array element-type, e.g. double-float.
There are some examples, and then the code.
I have also written some plotting functions (not included) so I can produce 1d plots, image
plots, surfaces, contours and filled contours. They work by producing, and then running, a
gnuplot script. You can do, e.g.
(plot (.sin (.* 2 th)) th)
Cheers, and thanks for the help.
David.
1d example showing how to create a 10 element array 'th' ranging between -pi and pi, and then
calculating sin(2*th)
[3]> (setf th (rseq (- pi) pi 10))
#(-3.1415926535897932385L0 -2.4434609527920614076L0 -1.7453292519943295769L0
-1.047197551196597746L0 -0.34906585039886591534L0 0.34906585039886591534L0
1.0471975511965977465L0 1.7453292519943295771L0 2.4434609527920614078L0
3.1415926535897932385L0)
[4]> (.sin (.* 2 th))
#(-1.0033115227368672049L-19 0.9848077530122080594L0 0.34202014332566873305L0
-0.8660254037844386469L0 -0.64278760968653932627L0 0.64278760968653932627L0
0.86602540378443864646L0 -0.34202014332566873345L0 -0.9848077530122080593L0
1.0033115227368672049L-19)
2d example that creates 2 10x10 arrays, in a similar fasion to Matlab's meshgrid, and then
calculates the
distance from the centre of the mesh to each point.
[8]> (setf x (rseq2 (- 1) 1 10))
#2A((-1 -1 -1 -1 -1 -1 -1 -1 -1 -1)
(-7/9 -7/9 -7/9 -7/9 -7/9 -7/9 -7/9 -7/9 -7/9 -7/9)
(-5/9 -5/9 -5/9 -5/9 -5/9 -5/9 -5/9 -5/9 -5/9 -5/9)
(-1/3 -1/3 -1/3 -1/3 -1/3 -1/3 -1/3 -1/3 -1/3 -1/3)
(-1/9 -1/9 -1/9 -1/9 -1/9 -1/9 -1/9 -1/9 -1/9 -1/9)
(1/9 1/9 1/9 1/9 1/9 1/9 1/9 1/9 1/9 1/9)
(1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3 1/3)
(5/9 5/9 5/9 5/9 5/9 5/9 5/9 5/9 5/9 5/9)
(7/9 7/9 7/9 7/9 7/9 7/9 7/9 7/9 7/9 7/9)
(1 1 1 1 1 1 1 1 1 1))
[9]> (setf y (transpose x))
#2A((-1 -7/9 -5/9 -1/3 -1/9 1/9 1/3 5/9 7/9 1)
(-1 -7/9 -5/9 -1/3 -1/9 1/9 1/3 5/9 7/9 1)
(-1 -7/9 -5/9 -1/3 -1/9 1/9 1/3 5/9 7/9 1)
(-1 -7/9 -5/9 -1/3 -1/9 1/9 1/3 5/9 7/9 1)
(-1 -7/9 -5/9 -1/3 -1/9 1/9 1/3 5/9 7/9 1)
(-1 -7/9 -5/9 -1/3 -1/9 1/9 1/3 5/9 7/9 1)
(-1 -7/9 -5/9 -1/3 -1/9 1/9 1/3 5/9 7/9 1)
(-1 -7/9 -5/9 -1/3 -1/9 1/9 1/3 5/9 7/9 1)
(-1 -7/9 -5/9 -1/3 -1/9 1/9 1/3 5/9 7/9 1)
(-1 -7/9 -5/9 -1/3 -1/9 1/9 1/3 5/9 7/9 1))
[10]> (.sqrt (.+ (.* x x) (.* y y)))
#2A((1.4142135
1.2668616
1.1439589
1.0540925...... [ Result cut because it's quite long...)
The code:
;; Calls function f for the ith element in the arrays in args
(defun f-one-element (f i args)
(let ((total (row-major-aref (car args) i)))
(dotimes (j (- (length args) 1))
(setf total (funcall f total (row-major-aref (nth (1+ j) args) i)))) total))
;; Calls function f for each element in the arrays in arg
;; and returns a new array
(defun array-f (f args)
(let ((r (make-array (array-dimensions (car args)))))
(dotimes (i (array-total-size (car args)))
(setf (row-major-aref r i) (f-one-element f i args) )) r))
;; Calls function f for each element in the arrays in arg
;; and the scaler in k
;; and returns a new array
(defun array-f-scaler (f k args)
(let ((r (make-array (array-dimensions (car args)))))
(dotimes (i (array-total-size (car args)))
(setf (row-major-aref r i) (funcall f (f-one-element f i args) k))) r))
(defun .f (f args)
(let ((arrays (remove-if-not (function arrayp) args))
(scalers (remove-if-not (function numberp) args)))
(if scalers
(array-f-scaler f (apply f (remove-if-not (function numberp) args)) arrays)
(array-f f arrays))))
;; Define add and multiply operators.
;; Can take any number of args.
(defun .* (&rest args)
(.f #'* args))
(defun .+ (&rest args)
(.f #'+ args))
;; Divide and minus take 1 or 2 args.
;; The order matters.
(defun .- (&rest args)
(let ((r (make-array (array-dimensions (car (remove-if-not (function arrayp) args))))))
(if (= (length args) 1)
(dotimes (i (array-total-size r))
(setf (row-major-aref r i) (- (row-major-aref (car args) i))))
(let ((a (car args))
(b (cadr args)))
(if (arrayp a)
(if (arrayp b)
(dotimes (i (array-total-size r))
(setf (row-major-aref r i)
(- (row-major-aref a i) (row-major-aref b i))))
(dotimes (i (array-total-size r))
(setf (row-major-aref r i)
(- (row-major-aref a i) b))))
(dotimes (i (array-total-size r))
(setf (row-major-aref r i)
(- a (row-major-aref b i))))))) r ))
;; Divide and minus take 1 or 2 args.
;; The order matters.
(defun ./ (&rest args)
(let ((r (make-array (array-dimensions (car (remove-if-not (function arrayp) args))))))
(if (= (length args) 1)
(dotimes (i (array-total-size r))
(setf (row-major-aref r i) (/ (row-major-aref (car args) i))))
(let ((a (car args))
(b (cadr args)))
(if (arrayp a)
(if (arrayp b)
(dotimes (i (array-total-size r))
(setf (row-major-aref r i)
(/ (row-major-aref a i) (row-major-aref b i))))
(dotimes (i (array-total-size r))
(setf (row-major-aref r i)
(/ (row-major-aref a i) b))))
(dotimes (i (array-total-size r))
(setf (row-major-aref r i)
(/ a (row-major-aref b i))))))) r ))
;; Array functions with 1 argument
(defun array-1f (f a)
(let ((r (make-array (array-dimensions a))))
(dotimes (i (array-total-size a))
(setf (row-major-aref r i) (funcall f (row-major-aref a i)))) r))
(defun .sin (a)
(array-1f #'sin a))
(defun .cos (a)
(array-1f #'cos a))
(defun .tan (a)
(array-1f #'tan a))
(defun .sqrt (a)
(array-1f #'sqrt a))
(defun .exp (a)
(array-1f #'exp a))
(defun .log (a) ; add &optional for base
(array-1f #'log a))
;;Vector generators
(defun iseq (start end)
(let ((r (make-array (- (1+ end) start))))
(dotimes (i (- (1+ end) start))
(setf (svref r i) i)) (.+ r start)))
(defun rseq (a b n)
(let ((r (make-array n))
(sp (/ (- b a) (1- n))))
(dotimes (i n)
(setf (svref r i) (+ (* i sp) a))) r))
(defun iseq2 (start end)
(let ((r (make-array (list (- (1+ end) start) (- (1+ end) start) ))))
(dotimes (i (- (1+ end) start))
(dotimes (j (- (1+ end) start))
(setf (aref r i j) i))) (.+ r start)))
(defun rseq2 (a b n)
(let ((r (make-array (list n n)))
(sp (/ (- b a) (1- n))))
(dotimes (i n)
(dotimes (j n)
(setf (aref r i j) (+ (* i sp) a)))) r))
(defun transpose (a)
(let ((r (make-array (list (array-dimension a 1) (array-dimension a 0)))))
(dotimes (i (array-dimension a 0))
(dotimes (j (array-dimension a 1))
(setf (aref r j i) (aref a i j))))r))
(defun rseqxy (a b n)
(values (rseq2 a b n) (transpose (rseq2 a b n))))
David Wallis <ยทยทยท@cpom.ucl.ac.uk> writes:
> ;; Calls function f for the ith element in the arrays in args
> (defun f-one-element (f i args)
> (let ((total (row-major-aref (car args) i)))
> (dotimes (j (- (length args) 1))
> (setf total (funcall f total (row-major-aref (nth (1+ j) args) i)))) total))
Please move the final "total))" to a separate line so that it
doesn't look like it's part of the loop. I'd also prefer longer
names for variables.
Instead of calling NTH repeatedly, you could use DOLIST:
(defun f-one-element (function index arrays)
(let ((total (row-major-aref (first arrays) index)))
(dolist (array (rest arrays))
(setf total (funcall function total
(row-major-aref array index))))
total))
On my machine, this change reduced the execution time by 10% for
#'* and three arrays.
The function could also be written with REDUCE, but that was
actually slower in the Lisp I'm using.