It's been a while since I wrote about my scientific computing library for lisp: map. Map is a great exercise for me -- in learning Lisp and learning algorithms and optimization. Since I'm between jobs and not yet in grad school, the time has been ripe to work on the library with some toy cases and try to make it as usable as possible by the time I get back in a position where I could use it for real work.

Because of my copious free time, a lot has been changing lately, including the addition of matrix slices, which were important, and such important features as matrix inversion. I'm also looking at adding parallelization at some point, as well as a matrix-defining macro in the style of the #v() one for vectors.

In any event, here's what's new:

Slices

I promised to write again when I had finished with matrix slices in map.  I'm sure they're not in their final state, but I have implemented them in a way that works. Now you can do things like:

(setf (mref m 1 #r(1 2) t) #2a((2.0d0 1.0d0 3.4d0) (3.1d0 3.2d2 4.2d1)))

Which already has a number of things going on in it.  This is setting a sub-matrix of m, a rank-3 matrix to a rank-2 matrix specified in the end of  the setf call.  Since we selected a single value for the first index of m, mref will not count that dimension as extant for the purposes of comparing with the shape of the matrix we're assigning.  This is really useful, because now things like multiplying matrices can be done like this:

(defmethod .* ((A simple-array) (B simple-array))
  (with-result (result (list (array-dimension A 0) (array-dimension B 1)))
    (do-matrix (result (i j))
      (setf (aref result i j) (dot (mref a i t) (mref b t j))))))

There are three different ways you can specify subscripts to mref:

  • Integers This will work like aref, and the dimension won't be included in the final output

  • Ranges This will select all entries in the matrix that have indices inside the specified range in this dimension, for example, (mref #(1 2 3) #r(0 1)) ==> #(1 2)
  • t This selects all elements in this dimension.  It's great for slicing rows/columns out of matrices a la (mref m t 0), which will take the first column out of m.

Matrix iteration

Since most of the functions in the library were operating on matrices and returning new ones, I thought some macros were in order.  I think it's the reflex of any lisper to write a macro when they see the same code in more than one place.

Of that urge were born (with-result) and (do-matrix).  (with-result) is simply a convenient way of defining a new matrix at the start of a function and automatically returning it.  (do-matrix), though, I rather like because it allows looping over matrices in a rather intuitive way.

(defun mandelbrot (matrix)
  (with-result (result (list 1001 1001) 'integer)
    (do-matrix (matrix (r i))
      (setf (aref result r i)
	    (loop
	       with c = (aref matrix r i)
	       with z = c
	       for i from 0 upto 80
	       if (> (abs z) 2)
	       return (1- i)
	       else
	       do (setf z (+ (expt z 2) c))
	       finally (return 80))))))

Here you can see (do-matrix) at work, testing whether each element of an array of complex numbers is in the Mandelbrot set.  It's definition is

(defmacro do-matrix (matrix &optional subscripts) ...)

which means you can then just write element-wise operations inside the loop with all of the indices figured out for you.

I'll end with the figure I made of the Mandelbrot set today on the [-2, 2] range in the complex plane. Along with the (mandelbrot) function above, you need the matrix to call it on,

(mute
  (setq cm
	(let ((center 500))
	  (with-result (result (list 1001 1001) '(complex double-float))
	    (do-matrix (matrix (r i))
	      (setf (aref matrix r i)
		    (complex (* (- r center) (/ 2.0d0 500))
			     (* (- i center) (/ 2.0d0 500)))))))))

and then you can run

(mute (setq r (mandelbrot cm)))
(image r)

.

{{< figure src="/wp-content/uploads/mandelbrot.png" title="Mandelbrot Set" alt="The Mandelbrot set from -2 to 2 in the complex plane." caption="The Mandelbrot set from -2 to 2 in the complex plane." >}}