07 September 2007
Quick Note on Continuations to comp.lang.scheme
15 August 2007
One More Reason to Love Module Systems and Hygiene
I just released an Automatic Differentiation PlaneT Package, and am now converting my numerical relativity code to use it. The package exports procedures like +
, 
, =
, etc, which are modified from the standard ones so it can compute derivatives of numerical procedures which use them. Of course, these modified procedures are slower than their standard cousins.
My numerical relativity code uses all sorts of looping constructs from the excellent SRFI42. These loops iterate over vectors, over ranges of numbers, lists, matrices, .... The iteration macros expand to code like (this is just an illustrationthe real code is more complicated than this)
(let loop ((vectorindex 0)) (if (>= vectorindex (vectorlength vector)) 'done (let ((vectorelt (vectorref vector vectorindex))) <body> (loop (+ vectorindex 1)))))
+
, >=
, etc. It is acceptable for these functions to be the addition and comparison from my PlaneT package, since that package preserves the meaning of the numerical functions operating on pure numbers. But, it would be nicer if the +
and >=
came from the mzscheme
language module. In this case, the compiler could inline the addition and comparison, and even if it didn't there wouldn't be the differentiation overhead.
No problem! PLT Scheme's module system takes care of this automatically: SRFI42 is defined in the mzscheme
language, so all references to +
and >=
in its macros refer to those bindings from mzscheme
. So, even if the module using the SRFI42 macros is written in a different language (like "derivlang.ss"
), the expansion of these macros uses the correct procedures.
This is a combined effect of of hygiene and the PLT module system. Neither alone would reproduce the correct behavior. For example, in schemes without a module system which handles macros properly, it would be necessary to redefine +
and >=
at the top level in the deriv.plt
package. Then, the hygiene of the macro expander would ensure that the +
and >=
in the macro expansion of SRFI42 iterators referred to the toplevel bindings (even if there were local bindings for these symbols at the macro usesite), but these would be the derivenabled +
and >=
, not the usual, optimizable mzscheme
ones.
In short: I know of no other module system that would handle this so elegantly (in a current scheme system, anywaythe proposed R6RS module system would do just as well for this). (In fact, I think only the module system of Chez scheme, immortalized in the portable syntaxcase expander by Dybvig and Waddell handles this at all. Unfortunately, it requires you to list the local variables which a macro can expand into manually along with the usual exports of the module, which is not nearly as convenient.) Hooray for PLT!
Automatic Differentiation PlaneT Package Released
08 August 2007
Persistency and Lazy Memoization in a Pairing Heap
I just released a pairing heap PlaneT package. Because the datastructure is both persistent and has amortized space bounds, it requires a neat trick to implement, which I will discuss below.
A pairing heap is an adaptive datastructure for a collection of objects which supports the following operations [edit: You can also read about an alternative implementation here in the Scheme Cookbook. This implementation doesn't achieve the amortized bounds on removemin
which mine does in the presence of persistence]:
Operation  Time Bound 
min  O(1) 
insert  O(1) 
merge  O(1) 
removemin  O(n) worst case, O(log(n)) amortized. 
The pairing heaps in my PlaneT package are persistent, which means that the above operations do not change the heap given to them, but rather construct a new heap which shares some structure with the old heap. For example, the (removemin h)
procedure returns a new heap which contains all the elements of h
except for the smallest element.
Note that the timebound on removemin in O(n) worstcase, and O(log(n)) amortized. This means that, while a single removemin operation can take time which is O(n), most will take much less. In fact, the longer removemin operations are so rare that, in a sequence of an infinite number of removemin operations, the average time per operation is O(log(n)). (This bound is only conjectured at the moment, but supported by a lot of experimental evidence. Because the pairingheap adapts its structure to the use pattern it observesa really neat feature, by the wayit is extremely difficult to rigorously prove bounds for it.)
It turns out to be difficult to ensure that amortized bounds still obtain with a persistent datastructure. The reason is illustrated in the following code snippet which creates a list of 1000 threads, each of which does something with a heap:
(let ((h (makeexpensiveremoveminheap))) (listec (:range i 1000) (thread (lambda () (let ((newh (removemin h))) (dosomethinginterestingwith newh))))))
removemin
here) can be "accounted" to a sequence of prior cheap operations (basically, makeexpensiveremoveminheap
must do a lot of inserts, each of which can "carry" a bit of the cost of the upcoming removemin
operation, so that you know that expensive removemin
s are very rare, and the cost averages out to O(log(n))). That style of accounting (often called the "Banker's Method") works if you know that only one removemin
will occur in the future of a given datastructure, as would be case if we were using ephemeral datastructuresin this case, the first removemin
would destroy the datastructure, so we could only perform it once. But we're not using ephemeral datastructures: each of the 1000 threads created above will reperform the expensive removemin
on the same h
, and we'll be in a pile of trouble, timeboundswise.
The solution to this problem is to use lazy evaluation with memoization. When constructing a heap, queue up the removemin
operation using delay
. When you need to remove the minimum element, use force
. Because delay
and force
memoize the result, you can guarantee that removemin
operations after the first will complete in O(1) time, which maintains the amortized bounds in the presence of persistence.
This is, of course, not my idea. It (and much more) are explained in the delightful thesis of Chris Okasaki, and also in his reportedlydelightful (I've never read it, but it's very wellregarded) book, Purely Functional Data Structures.
By the way, if you need a bunch of datastructures (more than my simple pairingheap), please check out the Galore PlaneT library by Carl Eastlund, Jens Axel Soegaard, and Stevie Strickland (either Version 2.2which contains more datastructures, but also more bugs and isn't as cleanor Version 3.6which is the most modern version). It's what I use when I need a serious finite map or set.
06 August 2007
I Just Voted "Yes" on R6RS
Yep, I did. I'm including my explanation below, because I think that R6RS has come under a lot of fire lately. (By the way, for those who say, "why doesn't the standard come with a small base, and then lots of libraries attached to it?", it does! <snark>Try reading it, why don't you.</snark>)
Even if you don't agree with my choice, don't forget to vote if you signed up (just one of many reminders you'll surely get until voting closes Sunday)!
My explanation:
I would gladly trade a report full of mistakesnot that I think this is such a report, but some dofor a library system which properly handles hygienic macros. Therefore, I am voting "Yes" on R6RS. (The exception specification is just icing on the cake, as far as I'm concerned.)
Some consolation for those who do think this report is full of mistakes:
 I expect we'll be seeing lots of "(r6rs base)"compliant implementations which leave off libraries they don't like, and

...in most matters it is more important that the applicable rule of law be settled than that it be settled right.
 Burnet v. Coronado Oil & Gas Co., 285 U.S. 393, 406408 (1932) (Justice Brandeis dissenting).
23 July 2007
A Question (and an Answer, sort of) About Scheme FFIs
Recently, I received the following email from Raoul Duke:
hi, I'm trying to learn up on what Scheme has the easiest FFI for talking to C (more ideally C++, but I assume nobody really does that). I've been reading up and it seems like PLT isn't bad. I came across your blog  might you have any opinions or recommendations about PLT vs. Chicken vs. Bigloo vs. Gambit vs. etc.? many thanks!Raoul generously allowed me to post his original mail here along with my response. I hope that it will be useful to others who have similar questions. I hope that I have my facts correct (I have some experience with this), but if I've made a mistake please correct me. Also, please feel free to post your opinions on this topic in the comments! My response:
Raoul,
I do have extensive experience with the FFI from PLT Scheme, and some experience with Chicken's, Bigloo's and Gambit's. Do you mind if I post this message and your original request to my blog (unless you request otherwise, I would leave out your email address, but include your name)? I think this is something that many other people would be interested in, too.
Essentially, the FFIs split into two types:
 PLT
 Access everything dynamically from Scheme (i.e. at runtime).
 Chicken, Bigloo, Gambit
 Write some Scheme code which, when *compiled*, generates C stubs that wrap a given C library.
(Chicken has a dynamic FFI, toosee the lazyffi egg herebut most people don't use it.)
Personally, I favor the PLT approach. The following is all you need to do, for example, to use the everpopular daxpy from the BLAS library (if you really want to do this, you should check out the pltlinalg PlaneT package):
(module demonstrateffi mzscheme (require (lib "foreign.ss")) (unsafe!) (define daxpy (getffiobj 'daxpy (ffilib "libblas") (_fun _int _double* _f64vector _int _f64vector _int > _void))))
Run this in DrScheme under the "module" language, and you've got daxpy. Of course, since you're doing all this binding at runtime, it's really easy to do all sorts of fancy stuff because you're in the Scheme world the whole time. See this paper for a discussion of the design that went into PLT's FFI. I cannot emphasize enough how great it is to be able to stay in the Scheme world when importing C functions.
The drawback: it doesn't do C++ at all. To bind to C++, you will have to write some code with extern "C" linkage to produce stub functions which you can bind using the FFI. That is possible to do, but it could be a pain, depending on how much you want to wrap. I think the PLT guys themselves did this with wxWindows in order to get their MrEd graphical framework; you might want to ask on the mailing list about this.
The alternate approach (used by Chicken, Bigloo, and Gambit) is not without advantages. Chicken, for one, can bind directly to C++ (at least to a large subset of C++); it imports C++ classes as tinyclos objects. And Chicken, in particular, provides a really nice C/C++ parser (see the easyffi egg here) that will automatically generate a lot of your bindings for you (which is something that PLT scheme will not do). Bigloo, also, has a tool called cigloo (though it's probably now named bglcigloo) which parses C header files and outputs the proper foreign declaration. I've had more trouble with it than with the Chicken toolsit gets confused easily. For both, I have found it's better to pipe or copy and paste only some pieces of the header into these tools so they don't get confused.
Gambit doesn't have a FFI parser (or at least I haven't found one for it); you just write the input and output types in a (clambda ...) special form.
Since all of Chicken, Bigloo, and Gambit require you to compile the FFI declarations before you can use them, you'll have to make sure that you have all the proper include and linking flags on the command line in order to make your library workthis can be a big drawback. In PLT, the only thing you need is to find the actual object file(s) which make up your library. You can do this using the extensive filesearching/existence testing, etc, code in PLT, and you can do it at runtime, so it's easy to keep searching for the file in lots of strange places, in case your users don't always put their libraries in /usr/local or whatever.
In short: I heartily recommend the PLT approach if you're willing to write some extern "C" stubs for your C++ code. If you're not, then I recommend Chicken, as it can (probably) handle the C++ natively. The only reason I see to use Bigloo or Gambit is if you require serious speed in the Scheme part of your code; in that case, use Bigloo if you are going to write C++ style code in Scheme, and Gambit if you need things like call/cc.
Hope this helps! Please feel free to ask more questions if any occur to you.
Will
17 June 2007
NBody Utilities Coming Real Soon Now (TM)
It's been a while since I posted here, but never fear: I've been busy. I've been working on an updated set of NBody utilities for another paper examining the efficiency of our algorithm (see astroph/0611416, to be published in ApJ 663 (2007) 1420) relative to driftkickdrift leapfrog in a treecode. I've got code up and running, and am beginning simulations to compare the two algorithms. The code is a mix of C (for anything O(N) and larger) and Scheme (to patch the various C algorithms together into actual simulations).
It looks, preliminarily, like our algorithm significantly outperforms leapfrog (i.e. treecomputed force errors dominate at larger timesteps with our algorithm) in this lowaccuracy regime, which is very exciting. Right now Ed has an idea for quantifying the errors introduced into this process in terms of the change in the velocity diffusion coefficient in the FokkerPlanck equation when the accuracy parameter of the treecode is varied, so I'm writing some code to run simulations and compute the velocity diffusion coefficient.
Once I've got the code into a workable state (i.e. once it easily runs all the simulations I've got set up for it), I'll write up some documentation and post the whole mess on PlaneT. NBody simulations for all!
I just read on this PLT Scheme mailing list thread about a nifty way to post syntaxhighlighted scheme code to your blog, so I have to try it out. Here's some code which runs a 10Kbody Plummer model up to core collapse (at which point, you'll have to use a code which models binary systems with small separations carefully, which this code does not).
(module uptocorecollapse mzscheme (require (planet "bhtree.ss" ("wmfarr" "nbody.plt" 1 0)) (planet "ics.ss" ("wmfarr" "nbody.plt" 1 0)) (planet "timestep.ss" ("wmfarr" "nbody.plt" 1 0)) (planet "evolve.ss" ("wmfarr" "nbody.plt" 1 0)) (planet "system.ss" ("wmfarr" "nbody.plt" 1 0)); (planet "gl3pt.ss" ("wmfarr" "nbody.plt" 1 0))) ;; The bodies. (define bs (makeplummer* 10000)) ;; Various parameters. (define T (* 20 (medianrelaxationtime bs))) ; Core collapse within 20 median relax times. (define eps 1e3) ; Force softening: don't resolve scales < 1e3 (define energybeta 0.1) (define advancebeta 1.0) (define outputinterval 1.0) ;; Total initial energy (define e (bhtreetotalenergy eps energybeta bs)) ;; Timestep function. (define frac 0.5) (define timestep (maxforcetimestep eps frac)) ;; Singlestep advancer. (define adv (lambda (bs dt) (gl3ptadvance eps advancebeta bs dt))) ;; Filter between stepsoutput roughly at outputinterval. (define lastoutput (makeparameter 0)) (define t (makeparameter 0)) (define (filter bs dt) (t (+ (t) dt)) (when (> (t) (+ (lastoutput) outputinterval)) (lastoutput (t)) (let ((enew (bhtreetotalenergy eps energybeta bs))) (printf (stringappend "At time ~a, energy is ~a, 10% lagrangian " "radius is ~a, energy error is ~a~%") (t) enew (lagrangianradius bs 0.1) (abs (/ ( e enew) e)))))) ;; Do it. ((evolve adv timestep filter) bs T))
29 April 2007
Pure Scheme Code for SO(3,1)
Jens Axel Soegaard wrote a comment on my last post asking for the Scheme code for comparison, so here it is. Note that the code provides a (moderately efficient, but conses too much) matrix inversion routine in pure Scheme (you want to look at the matrix module, for matrixinverse)some people may find that useful independently. Sorry it's so long, but there's a lot of background stuff to get through (it's part of a larger code base, too, but I think this is enough to be selfcontained). Enjoy!
Here's the code for the generics that I use at the end of the matrix code:
# base.ss: Generic operations and base language. Copyright (C) 2007 Will M. FarrThis program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 USA. # (module base (lib "swindle.ss" "swindle") (provide (allfromexcept (lib "swindle.ss" "swindle") +  * / add negate) (rename + mz:+) (rename  mz:) (rename * mz:*) (rename / mz:/) (rename my+ +) (rename my ) (rename my* *) (rename my/ /) (rename myadd add) sub mul div invert (rename mynegate negate)) (define my+ (caselambda (() 0) ((x) x) ((x y) (myadd x y)) ((x y . ys) (apply my+ (myadd x y) ys)))) (define my (caselambda ((x) (mynegate x)) ((x y) (sub x y)) ((x . ys) (sub x (apply my+ ys))))) (define my* (caselambda (() 1) ((x) x) ((x y) (mul x y)) ((x y . ys) (apply my* (mul x y) ys)))) (define my/ (caselambda ((x) (invert x)) ((x y) (div x y)) ((x . ys) (div x (apply my* ys))))) (defgeneric (myadd x y)) (defmethod (myadd (x ) (y )) (+ x y)) (defgeneric (mynegate x)) (defmethod (mynegate (x )) ( x)) (defgeneric (sub x y)) (defmethod (sub (x ) (y )) ( x y)) (defgeneric (mul x y)) (defmethod (mul (x ) (y )) (* x y)) (defgeneric (invert x)) (defmethod (invert (x )) (/ x)) (defgeneric (div x y)) (defmethod (div (x ) (y )) (/ x y)))
Here's the matrix module (before we get to the Lorentz group):
# matrix.ss: Basic matrix utilities. Copyright (C) Will M. FarrThis program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 021101301, USA. # (module matrix mzscheme (require (lib "42.ss" "srfi") (lib "contract.ss") (only "base.ss" structtype>class add sub mul div invert negate defmethod )) (define mz:+ +) (define mz: ) (define mz:* *) (define mz:/ /) (definestruct matrix (rows cols elts) #f) (define (structtype>class struct:matrix)) (provide struct:matrix :matrix matrixec *small* ) (define (validmatrixrow/c m) (and/c naturalnumber/c (between/c 0 (mz: (matrixrows m) 1)))) (define (validmatrixcol/c m) (and/c naturalnumber/c (between/c 0 (mz: (matrixcols m) 1)))) (define (samematrixdimensions/c m) (flatnamedcontract "same matrix dimensions" (let ((mc (matrixcols m)) (mr (matrixrows m))) (lambda (m2) (and (= (matrixcols m2) mc) (= (matrixrows m2) mr)))))) (define squarematrix/c (flatnamedcontract "square matrix" (lambda (m) (= (matrixrows m) (matrixcols m))))) (provide/contract (rename mymakematrix makematrix (> naturalnumber/c naturalnumber/c any/c matrix?)) (rename mymatrix matrix (>r ((r naturalnumber/c) (c naturalnumber/c)) elts (and/c list? (lambda (els) (= (length els) (mz:* r c)))) matrix?)) (matrix? (> any/c boolean?)) (matrixref (>r ((m matrix?) (i (validmatrixrow/c m)) (j (validmatrixcol/c m))) any)) (matrixset! (>r ((m matrix?) (i (validmatrixrow/c m)) (j (validmatrixcol/c m)) (x any/c)) any)) (matrixmul (>r ((m1 matrix?) (m2 (and/c matrix? (flatnamedcontract "compatible matrix dimensions" (lambda (m2) (= (matrixcols m1) (matrixrows m2))))))) matrix?)) (matrixtranspose (> matrix? matrix?)) (matrixvectormul (>r ((m matrix?) (v (and/c vector? (flatnamedcontract "compatible vector dimension" (lambda (v) (= (vectorlength v) (matrixcols m))))))) vector?)) (vectormatrixmul (>r ((v vector?) (m (and/c matrix? (flatnamedcontract "compatible vector dimension" (lambda (m) (= (vectorlength v) (matrixrows m))))))) vector?)) (matrixadd (>r ((m1 matrix?) (m2 (and/c matrix? (samematrixdimensions/c m1)))) matrix?)) (matrixsub (>r ((m1 matrix?) (m2 (and/c matrix? (samematrixdimensions/c m1)))) matrix?)) (matrixscale (>r ((s (or/c number? matrix?)) (m (if (number? s) matrix? number?))) matrix?)) (matrixnorm (> (and/c integer? (>/c 0)) (> matrix? number?))) (matrixinfinitynorm (> matrix? number?)) (matrixinverse (> squarematrix/c squarematrix/c)) (matrixidentity (> naturalnumber/c matrix?)) (matrixcopy (> matrix? matrix?))) (define (mymakematrix rows cols elt) (makematrix rows cols (makevector (mz:* rows cols) elt))) (define (matrixref m i j) (vectorref (matrixelts m) (mz:+ (mz:* i (matrixcols m)) j))) (define (matrixset! m i j x) (vectorset! (matrixelts m) (mz:+ (mz:* i (matrixcols m)) j) x)) (define (mymatrix rows cols . elts) (makematrix rows cols (list>vector elts))) (definesyntax matrixec (syntaxrules () ((matrixec rrows ccols etc ...) (let ((rows rrows) (cols ccols)) (makematrix rows cols (vectoroflengthec (mz:* rows cols) etc ...)))))) (definesyntax :matrix (syntaxrules (index) ((:matrix cc var arg) (:vector cc var (matrixelts arg))) ((:matrix cc var (index i j) arg) (:do cc (let ((m arg) (rows #f) (cols #f)) (set! rows (matrixrows m)) (set! cols (matrixcols m))) ((i 0) (j 0)) (< i rows) (let ((i+1 (mz:+ i 1)) (j+1 (mz:+ j 1)) (wrapping? #f)) (set! wrapping? (>= j+1 cols))) #t ((if wrapping? i+1 i) (if wrapping? 0 j+1)))))) (define (matrixmul m1 m2) (let ((rows (matrixrows m1)) (cols (matrixcols m2)) (n (matrixcols m1))) (matrixec rows cols (:range i rows) (:range j cols) (sumec (:range k n) (* (matrixref m1 i k) (matrixref m2 k j)))))) (define (matrixtranspose m) (matrixec (matrixcols m) (matrixrows m) (:range i (matrixcols m)) (:range j (matrixrows m)) (matrixref m j i))) (define (matrixvectormul m v) (vectoroflengthec (matrixrows m) (:range i (matrixrows m)) (sumec (:range j (matrixcols m)) (* (matrixref m i j) (vectorref v j))))) (define (vectormatrixmul v m) (vectoroflengthec (matrixcols m) (:range j (matrixcols m)) (sumec (:range i (matrixrows m)) (* (vectorref v i) (matrixref m i j))))) (define (matrixadd m1 m2) (matrixec (matrixrows m1) (matrixcols m1) (:parallel (:matrix x1 m1) (:matrix x2 m2)) (+ x1 x2))) (define (matrixsub m1 m2) (matrixec (matrixrows m1) (matrixcols m1) (:parallel (:matrix x1 m1) (:matrix x2 m2)) ( x1 x2))) (define (matrixscale m s) (when (number? m) (let ((temp m)) (set! m s) (set! s temp))) (matrixec (matrixrows m) (matrixcols m) (:matrix x m) (* x s))) (define ((matrixnorm n) m) (expt (sumec (:matrix x m) (expt (abs x) n)) (/ n))) (define (matrixinfinitynorm m) (maxec (:matrix x m) (abs x))) (define (scalematrixrow! m i s) (doec (:range j (matrixcols m)) (matrixset! m i j (* (matrixref m i j) s)))) (define swapmatrixrows! (caselambda ((m i j pvec) (let ((oldi (vectorref pvec i))) (vectorset! pvec i (vectorref pvec j)) (vectorset! pvec j oldi)) (swapmatrixrows! m i j)) ((m i j) (doec (:range k (matrixcols m)) (let ((oldik (matrixref m i k))) (matrixset! m i k (matrixref m j k)) (matrixset! m j k oldik)))))) (define (addscaledmatrixrowtorow! m s i j) (doec (:range k (matrixcols m)) (matrixset! m j k (+ (matrixref m j k) (* s (matrixref m i k)))))) (define (pivotrow m j) (fold3ec 'nomax (:range i j (matrixrows m)) i (lambda (x) x) (lambda (i iold) (if (> (abs (matrixref m i j)) (abs (matrixref m iold j))) i iold)))) (define (matrixpref m p i j) (matrixref m (vectorref p i) j)) (define (matrixidentity n) (matrixec n n (:range i n) (:range j n) (if (= i j) 1 0))) (define (matrixcopy m) (matrixec (matrixrows m) (matrixcols m) (:matrix x m) x)) (define *small* (makeparameter 1.4901161193847656e08 (lambda (sm) (if (not (>= sm 0)) (error '*small* "Cannot set to negative number: ~a" sm) sm)))) (define (matrixinverse min) (let ((rows (matrixrows min)) (cols (matrixcols min))) (let ((m (matrixcopy min)) (id (matrixidentity rows)) (pvec (vectoroflengthec rows (:range i rows) i))) (doec (:range i rows) (begin (let ((pi (pivotrow m i))) (swapmatrixrows! m i pi pvec) (swapmatrixrows! id i pi)) (let ((pelt (matrixref m i i))) (when (< (abs pelt) (*small*)) (raisemismatcherror 'matrixinverse "Singular or nearly singular (by *small* parameter) matrix " (cons pelt min))) (scalematrixrow! m i (/ pelt)) (scalematrixrow! id i (/ pelt)) (doec (:range j rows) (not (= j i)) (let ((s ( (matrixref m j i)))) (addscaledmatrixrowtorow! m s i j) (addscaledmatrixrowtorow! id s i j)))))) id))) (define (relativesize x norm) (abs (/ x norm))) (defmethod (add (m1 ) (m2 )) (matrixadd m1 m2)) (defmethod (sub (m1 ) (m2 )) (matrixsub m1 m2)) (defmethod (mul (m ) (s )) (matrixscale m s)) (defmethod (mul (s ) (m )) (matrixscale m s)) (defmethod (mul (m ) (v )) (matrixvectormul m v)) (defmethod (mul (v ) (m )) (vectormatrixmul v m)) (defmethod (mul (m1 ) (m2 )) (matrixmul m1 m2)) (defmethod (negate (m )) (matrixscale m 1)) (defmethod (invert (m )) (matrixinverse m)) (defmethod (div (v ) (m )) (matrixvectormul (matrixinverse m) v)) (defmethod (div (m ) (s )) (matrixscale m (mz:/ s))))
And here's (finally), the code for SO(3,1)see the documentation for SO31.plt for an explanation of how we parameterize the group if this seems a little opaque :).
# lorentztransform.ss: SO(3,1) matrix group. Copyright (C) 2007 Will M. Farr This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 021101301, USA. # (module lorentztransform mzscheme (require "matrix.ss" (lib "42.ss" "srfi") (lib "contract.ss") (lib "math.ss") "simplex.ss" "simplexmaputils.ss" (only "base.ss" structtype>class defmethodmul div invert) (allexcept (planet "table.ss" ("soegaard" "galore.plt" 3)) map ormap andmap foreach values equal? union)) (definestruct lorentztransformfield (table) #f) (definestruct lorentztransformelement (params inverseparams matrix inversematrix) #f) (define (structtype>class struct:lorentztransformfield)) (define (structtype>class struct:lorentztransformelement)) (provide eta ltparams/c struct:lorentztransformfield struct:lorentztransformelement ) (define ltparams/c (and/c vector? (flatnamedcontract "6element vector" (lambda (v) (= (vectorlength v) 6))))) (define zerosimplex/c (and/c simplex? (flatnamedcontract "simplexdimension = 0" (lambda (s) (= (simplexdimension s) 0))))) (define ltfalist/c (listof (cons/c zerosimplex/c (or/c matrix? ltparams/c lorentztransformelement?)))) (provide/contract (matrixetatrace (> matrix? number?)) (lorentzmatrix (> number? number? number? number? number? number? matrix?)) (makelorentzmatrix (> ltparams/c matrix?)) (lorentzmatrix>params (> matrix? ltparams/c)) (Dlorentzmatrix (> ltparams/c (vectorof matrix?))) (lorentztransformfield? (> any/c boolean?)) (lorentztransformelement? (> any/c boolean?)) (params>lorentztransformelement (> ltparams/c lorentztransformelement?)) (matrix>lorentztransformelement (> matrix? lorentztransformelement?)) (lorentztransformelementparams (> lorentztransformelement? ltparams/c)) (lorentztransformelementinverseparams (> lorentztransformelement? ltparams/c)) (lorentztransformelementmatrix (> lorentztransformelement? matrix?)) (lorentztransformelementinversematrix (> lorentztransformelement? matrix?)) (lorentztransformelementcompose (> lorentztransformelement? lorentztransformelement? lorentztransformelement?)) (lorentztransformelementinvert (> lorentztransformelement? lorentztransformelement?)) (lorentztransformelementtransformvector (> lorentztransformelement? vector? vector?)) (rename mymakelorentztransformfield makelorentztransformfield (> ltfalist/c lorentztransformfield?)) (rename mylorentztransformfield lorentztransformfield (>* () ltfalist/c (lorentztransformfield?))) (lorentztransformfieldelement (> lorentztransformfield? zerosimplex/c lorentztransformelement?)) (lorentztransformfieldcompose (> lorentztransformfield? lorentztransformfield? lorentztransformfield?)) (lorentztransformfieldfold (> (> zerosimplex/c lorentztransformelement? any/c any) any/c lorentztransformfield? any))) (define (matrixetatrace M) (+ ( (matrixref M 0 0)) (matrixref M 1 1) (matrixref M 2 2) (matrixref M 3 3))) (define Rxgen (matrix 4 4 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0)) (define Rygen (matrix 4 4 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0)) (define Rzgen (matrix 4 4 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0)) (define Bxgen (matrix 4 4 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0)) (define Bygen (matrix 4 4 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0)) (define Bzgen (matrix 4 4 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0)) (define eta (matrix 4 4 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1)) (define (Rx x) (matrix 4 4 1 0 0 0 0 1 0 0 0 0 (cos x) ( (sin x)) 0 0 (sin x) (cos x))) (define (Ry y) (matrix 4 4 1 0 0 0 0 (cos y) 0 ( (sin y)) 0 0 1 0 0 (sin y) 0 (cos y))) (define (Rz z) (matrix 4 4 1 0 0 0 0 (cos z) ( (sin z)) 0 0 (sin z) (cos z) 0 0 0 0 1)) (define (Bx x) (matrix 4 4 (cosh x) ( (sinh x)) 0 0 ( (sinh x)) (cosh x) 0 0 0 0 1 0 0 0 0 1)) (define (By y) (matrix 4 4 (cosh y) 0 ( (sinh y)) 0 0 1 0 0 ( (sinh y)) 0 (cosh y) 0 0 0 0 1)) (define (Bz z) (matrix 4 4 (cosh z) 0 0 ( (sinh z)) 0 1 0 0 0 0 1 0 ( (sinh z)) 0 0 (cosh z))) (define gens (list Bxgen Bygen Bzgen Rxgen Rygen Rzgen)) (define matfuns (list Bx By Bz Rx Ry Rz)) (define rgens (reverse gens)) (define rmatfuns (reverse matfuns)) (define (makelorentzmatrix paramv) (let ((! vectorref) (* matrixmul) (p paramv)) (* (Bx (! p 0)) (* (By (! p 1)) (* (Bz (! p 2)) (* (Rx (! p 3)) (* (Ry (! p 4)) (Rz (! p 5))))))))) (define (lorentzmatrix vx vy vz x y z) (makelorentzmatrix (vector vx vy vz x y z))) (define (Dlorentzmatrix params) (let* ((lparams (vector>list params)) (mats (map (lambda (f x) (f x)) matfuns lparams)) (rmats (reverse mats))) (vectoroflengthec 6 (:range i 6) (foldec (matrixidentity 4) (:parallel (:list M rmats) (:list g rgens) (:range j 5 1 1)) (if (= i j) (matrixmul M g) M) matrixmul)))) (define (asinh z) (log (+ z (sqrt (+ (sqr z) 1))))) (define (acosh z) (log (+ z (* (sqrt ( z 1)) (sqrt (+ z 1)))))) (define (lorentzmatrix>params M) (let* ((sinhvz ( (matrixref M 3 0))) (vz (asinh sinhvz)) (coshvz (cosh vz)) (sinhvy ( (/ (matrixref M 2 0) coshvz))) (vy (asinh sinhvy)) (coshvy (cosh vy)) (sinhvx ( (/ (matrixref M 1 0) (* coshvy coshvz)))) (vx (asinh sinhvx))) (let ((newM (matrixmul (Bz ( vz)) (matrixmul (By ( vy)) (matrixmul (Bx ( vx)) M))))) (let* ((siny ( (matrixref newM 1 3))) (y (asin siny)) (cosy (cos y)) (sinz ( (/ (matrixref newM 1 2) cosy))) (z (asin sinz)) (cosz (cos z)) (cosx (/ (matrixref newM 3 3) cosy)) (sinx (/ (+ (matrixref newM 3 2) (* siny sinz cosx)) cosz)) (x (asin sinx))) (vector vx vy vz x y z))))) (define (mylorentztransformfield . maps) (mymakelorentztransformfield maps)) (define (params>lorentztransformelement p) (let* ((M (makelorentzmatrix p)) (Minv (matrixinverse M)) (pinv (lorentzmatrix>params Minv))) (makelorentztransformelement p pinv M Minv))) (define (matrix>lorentztransformelement M) (let* ((Minv (matrixinverse M)) (p (lorentzmatrix>params M)) (pinv (lorentzmatrix>params Minv))) (makelorentztransformelement p pinv M Minv))) (define (mymakelorentztransformfield ltfalist) (let ((eltalist (map (lambda (sporM) (let ((porM (cdr sporM))) (cons (car sporM) (cond ((vector? porM) (params>lorentztransformelement porM)) ((matrix? porM) (matrix>lorentztransformelement porM)) (else porM))))) ltfalist))) (makelorentztransformfield (alist>ordered simplexcompare eltalist)))) (define (lorentztransformfieldelement ltf s) (lookup s (lorentztransformfieldtable ltf) (lambda () (makelorentztransformelement (makevector 6 0) (makevector 6 0) (matrixidentity 4) (matrixidentity 4))) values)) (define (lorentztransformelementcompose elt1 elt2) (let* ((M1 (lorentztransformelementmatrix elt1)) (M2 (lorentztransformelementmatrix elt2)) (M (matrixmul M1 M2)) (p (lorentzmatrix>params M)) (Minv (matrixinverse M)) (pinv (lorentzmatrix>params Minv))) (makelorentztransformelement p pinv M Minv))) (defmethod (mul (e1 ) (e2 )) (lorentztransformelementcompose e1 e2)) (define (lorentztransformelementtransformvector elt v) (matrixvectormul (lorentztransformelementmatrix elt) v)) (defmethod (mul (e ) (v )) (lorentztransformelementtransformvector e v)) (define (lorentztransformelementinvert elt) (makelorentztransformelement (lorentztransformelementinverseparams elt) (lorentztransformelementparams elt) (lorentztransformelementinversematrix elt) (lorentztransformelementmatrix elt))) (defmethod (invert (e )) (lorentztransformelementinvert e)) (defmethod (div (v ) (e )) (matrixvectormul (lorentztransformelementinversematrix e) v)) (define lorentztransformfieldcompose (let ((tablecompose (simplexmapliftcombiners (lambda () (makeordered simplexcompare)) insert (lambda (sl el sr er accu) (insert sl (lorentztransformelementcompose el er) accu)) insert))) (lambda (ltf1 ltf2) (makelorentztransformfield (tablecompose (lorentztransformfieldtable ltf1) (lorentztransformfieldtable ltf2)))))) (defmethod (mul (ltf1 ) (ltf2 )) (lorentztransformfieldcompose ltf1 ltf2)) (define (lorentztransformfieldfold f start ltf) (fold f start (lorentztransformfieldtable ltf))))
28 April 2007
Why Use the FFI for Simple Numerical Operations in MzScheme?
Brad Lucier (of Numerical Partial Differential Equations in Scheme fame) posted a comment on my last post announcing the code for numerical manipulation of Lorentz group elements which reads:
I'm surprised that you decided to use the C FFI to manipulate 4x4 matrices; they are so small that it would seem better just to use mzscheme's numerical facilities which are getting better all the time.
This is an interesting question to answer, so I'm going to devote this post to answering it instead of leaving a short response in the comments section.
Before writing this code (which uses some C to do the manipulations), I wrote some code which used purely Scheme to perform the following sequence of operations, given a set of parameters for an SO(3,1) element:
 Compute the matrix associated with the parameters, M.
 Invert that matrix to produce another matrix Minv.
 Compute the parameters associated with Minv, pinv.
Code Style  Total CPU Time  GC Time 
Pure Scheme  51 sec  27 sec 
Mixed Scheme/C  4 sec  2.8 sec 
To philosophize: the role I see for Scheme in my computational world is as an excellent driver to glue a bunch of simple C or Fortran computational elements together. (This isn't my ideahave a look at this for a similar approach with Python.) PLT Scheme is my implementation of choice for this because of its excellent FFI, it's extensive libraries (particularly its support of software contracts), and it's relatively easy C functions for memory allocation and Scheme object manipulation. The contracts, in particular, are a big help, because they let me avoid doing any argument checking on the C side (where I want to minimize coding). Memory allocation is nice, because with a combination of mzc xform
and scheme_malloc[_atomic]
, I never again have to mess with malloc
and free
in the C that I write. I've written this code with an eye toward including it in my thesis work, which will also be written in this style, so it seemed natural to do the manipulations in C.
This is, of course, my view. People like Brad, Marc, Manuel, and Jeff have a different vision, and have done tremendous work in optimizing Scheme performance on numeric codes. I've tried that approach, and always find myself back with PLT's scheme and C.
SO(3,1) (Lorentz Group) Manipulations from PLT Scheme
I just released a PlaneT package which manipulates elements of the Lorentz group (here). Get it using
(require (planet "SO31.ss" ("wmfarr" "SO31.plt" 1 0)))At the moment, it only runs on Mac OS 10.3 or 10.4, ppc, because it uses some C code I wrote (for efficiency), and I haven't had the patience to understand all the hoops I need to jump through to get PLT Scheme to compile this to a library on every system.
If someone on Windows or Linux (or an x86 mac) really wants to use the package, it's pretty easy to compile manually (after downloading and installing from PlaneT): just compile the SO31.c file in the SO31c subdirectory of the package to a dynamic library in the same directory named "libSO31.
I hope people find all this linear algebra stuff useful (it certainly will be for my thesis code).
26 April 2007
Linear Algebra in PLT Scheme
Last night, I posted a library which uses BLAS and LAPACK from within PLT Scheme to perform simple numerical linear algebra.
In doing so, I've learned a lot about the FFI for PLT Scheme. The conclusion: it's great! In particular, the "custom _fun types", discussed here are tremendously helpful! Particularly when interfacing to Fortran libraries, where every argument to a function must be a pointer. For example, here's how I make a call to dgesv
from LAPACK to solve linear systems:
(define matrixsolvemany (getffiobj 'dgesv_ *lapack* (_fun (m b) :: ((_ptr i _int) = (matrixrows m)) ((_ptr i _int) = (matrixcols b)) (_matrix = (matrixcopy m)) ((_ptr i _int) = (matrixrows m)) (_u32vector o (matrixrows m)) (x : _matrix = (matrixcopy b)) ((_ptr i _int) = (matrixrows b)) (_ptr o _int) > _void > x)))all the
(_ptr o ...)
calls allocate the appropriate memory blocks and pass pointers to them for dgesv
, and then handle the translation back to the scheme world automatically. I didn't have to write any C code to construct the library.
Up next: C code and PLT FFI for fast manipulation of SO(3,1) group elements (i.e. Lorentz transformations). Stay tuned.
09 April 2007
Pearls Before Breakfast
03 March 2007
Eager Comprehensions for Arrays
# arrayec.ss: SRFI42 comprehensions for SRFI25 arrays. Copyright (C) 2007 Will M. FarrThis library is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation; either version 2.1 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with this library; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 USA # (module arrayec mzscheme (require (lib "42.ss" "srfi") (lib "25.ss" "srfi") (only (lib "1.ss" "srfi") fold)) (provide (allfrom (lib "42.ss" "srfi")) (allfrom (lib "25.ss" "srfi")) arrayec :array) (definesyntax arrayec (syntaxrules () ((arrayec shp expr expr2 ... finalexpr) (let ((a (makearray shp))) (let ((arank (arrayrank a))) (let ((sizes (listec (:range i arank) ( (arrayend a i) (arraystart a i))))) (let ((indexmap (lambda (i) (apply values (car (fold (lambda (size indicesandi) (let ((indices (car indicesandi)) (i (cdr indicesandi))) (cons (cons (modulo i size) indices) (quotient i size)))) (cons '() i) (reverse sizes))))))) (let ((shareda (sharearray a (shape 0 (productec (:list size sizes) size)) indexmap))) (let ((i 0)) (doec expr expr2 ... (begin (arrayset! shareda i finalexpr) (set! i (+ i 1)))) a))))))))) (define (makeindexgenerator arr) (let ((r (arrayrank arr))) (let ((r1 ( r 1)) (lbs (vectoroflengthec r (:range i r) (arraystart arr i))) (ubs (vectoroflengthec r (:range i r) (arrayend arr i)))) (let ((idxv (vectoroflengthec r (:vector lb lbs) lb))) (vectorset! idxv r1 ( (vectorref idxv r1) 1)) (lambda () (vectorset! idxv r1 (+ (vectorref idxv r1) 1)) (let loop ((i r1)) (cond ((= i 0) (if (>= (vectorref idxv 0) (vectorref ubs 0)) 'done idxv)) ((= (vectorref idxv i) (vectorref ubs i)) (vectorset! idxv i (vectorref lbs i)) (let ((i1 ( i 1))) (vectorset! idxv i1 (+ (vectorref idxv i1) 1)) (loop i1))) (else idxv)))))))) (definesyntax :array (syntaxrules (index) ((:array cc x (index k0 ...) arrexpr) (:do cc (let ((arr arrexpr) (gen #f)) (set! gen (makeindexgenerator arr))) ((idxv (gen))) (not (eq? idxv 'done)) (let ((i 0) (x (arrayref arr idxv)) (k0 #f) ...) (begin (set! k0 (vectorref idxv i)) (set! i (+ i 1))) ...) #t ((gen)))) ((:array cc x arrexpr) (:array cc x (index) arrexpr)))))
18 February 2007
BarnesHut Tree Code Submitted to PlaneT
This is just another step in my code cleanup in preparation for comparing our algorithm to the standard leapfrog algorithm of cosmological codes.
NBody Simulation Initial Conditions
The paper is undergoing its final revisions, and the referee has indicated that these will make it acceptable, so I think we're in the home stretch of getting it published. That means that I'm consolidating my code a bit in preparation for testing our algorithm on BarnesHutstyle cosmological simulations. The PlaneT package is a result of that consolidation.
My current plan is to gin up a prototype BarnesHut tree code in PLT Scheme (shouldn't take longit's really pretty easy to code). Once I've got that working, I'll start converting the important bits to C that I can link into PLT Scheme. This is for speed. The IC code currently takes about 2 minutes on my machine to generate a millionbody system, while similar C code I've used takes as short as 10 seconds. There will likely be a similar factor for the integrator code.
We'll be wanting to test our integrator on large systems (i.e. around a million bodies), so coding the bottlenecks in C will be necessary. I'll be using as little of it as possible, however, because scheme is so much more pleasant.
I'm looking at this future work as a grand experiment to see whether it's practical to develop physics simulation code this way. I know that other people have done it before, but I want to see how it works for me.