Scheming: 2007

07 September 2007

Quick Note on Continuations to comp.lang.scheme

I just posted a quick note on continuations to c.l.s. I'm afraid it's pretty basic, and doesn't really do the topic justice, but maybe the original poster will find it helpful. I'm noting it here because no ever seems to explain continuations the way I did in my note (at least, I haven't seen any explanations like this), which uses "return points" to illustrate what continuations do. I don't fully understand continuations (in particular, I have a lot of trouble with the delimited continuations that PLT Scheme has recently acquired), but thinking of them the way I talk about in the c.l.s. post is the way I manage to muddle through using them. How do others think about call/cc? (i.e. what images are going through your head while you figure out how to do what you want with call/cc in code, not how do you formalize call/cc, or how you explain it to someone else, or ....)

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 SRFI-42. These loops iterate over vectors, over ranges of numbers, lists, matrices, .... The iteration macros expand to code like (this is just an illustration---the real code is more complicated than this)

(let loop ((vector-index 0))
  (if (>= vector-index (vector-length vector))
      'done
      (let ((vector-elt (vector-ref vector vector-index)))
        <body>
        (loop (+ vector-index 1)))))

You'll note that this code uses +, >=, 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: SRFI-42 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 SRFI-42 macros is written in a different language (like "deriv-lang.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 re-define + 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 SRFI-42 iterators referred to the top-level bindings (even if there were local bindings for these symbols at the macro use-site), but these would be the deriv-enabled + 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, anyway---the 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 syntax-case 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

I just released a PlaneT Package which, when used as a PLT Scheme module's initial language import, allows to compute derivatives of functions written in the module exactly. (I've written about the technique here and in the posts which follow.) I'm going to use this code to compute some derivatives of gravitational actions in code which would be extremely complicated to do by hand. Enjoy!

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 remove-min which mine does in the presence of persistence]:

Operation	Time Bound
min	O(1)
insert	O(1)
merge	O(1)
remove-min	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 (remove-min h) procedure returns a new heap which contains all the elements of h except for the smallest element.

Note that the time-bound on remove-min in O(n) worst-case, and O(log(n)) amortized. This means that, while a single remove-min operation can take time which is O(n), most will take much less. In fact, the longer remove-min operations are so rare that, in a sequence of an infinite number of remove-min 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 pairing-heap adapts its structure to the use pattern it observes---a really neat feature, by the way---it 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 (make-expensive-remove-min-heap)))
  (list-ec (:range i 1000)
    (thread
     (lambda ()
       (let ((new-h (remove-min h)))
         (do-something-interesting-with new-h))))))

A persistent datastructure can have multiple futures! But, when proving amortized bounds, one typically exploits that the cost of an expensive operation (remove-min here) can be "accounted" to a sequence of prior cheap operations (basically, make-expensive-remove-min-heap must do a lot of inserts, each of which can "carry" a bit of the cost of the upcoming remove-min operation, so that you know that expensive remove-mins 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 remove-min will occur in the future of a given datastructure, as would be case if we were using ephemeral datastructures---in this case, the first remove-min 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 re-perform the expensive remove-min on the same h, and we'll be in a pile of trouble, time-bounds-wise.

The solution to this problem is to use lazy evaluation with memoization. When constructing a heap, queue up the remove-min operation using delay. When you need to remove the minimum element, use force. Because delay and force memoize the result, you can guarantee that remove-min 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 reportedly-delightful (I've never read it, but it's very well-regarded) book, Purely Functional Data Structures.

By the way, if you need a bunch of datastructures (more than my simple pairing-heap), please check out the Galore PlaneT library by Carl Eastlund, Jens Axel Soegaard, and Stevie Strickland (either Version 2.2---which contains more datastructures~~, but also more bugs and isn't as clean~~---or Version 3.6---which 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 mistakes---not that I think this is such a report, but some do---for 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, 406--408 (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, too---see the lazy-ffi egg here---but 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 ever-popular daxpy from the BLAS library (if you really want to do this, you should check out the plt-linalg PlaneT package):

(module demonstrate-ffi mzscheme
  (require (lib "foreign.ss"))

  (unsafe!)

  (define daxpy
    (get-ffi-obj 'daxpy (ffi-lib "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 tools---it 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 (c-lambda ...) 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 work---this 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 file-searching/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

N-Body 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 N-Body utilities for another paper examining the efficiency of our algorithm (see astro-ph/0611416, to be published in ApJ 663 (2007) 1420) relative to drift-kick-drift 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. tree-computed force errors dominate at larger timesteps with our algorithm) in this low-accuracy 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 Fokker-Planck 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. N-Body simulations for all!

I just read on this PLT Scheme mailing list thread about a nifty way to post syntax-highlighted scheme code to your blog, so I have to try it out. Here's some code which runs a 10K-body 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 up-to-core-collapse mzscheme
  (require (planet "bh-tree.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 (make-plummer* 10000))
  
  ;; Various parameters.
  (define T (* 20 (median-relaxation-time bs))) ; Core collapse within 20 median relax times.
  (define eps 1e-3) ; Force softening: don't resolve scales < 1e-3
  (define energy-beta 0.1)
  (define advance-beta 1.0)
  (define output-interval 1.0)
  
  ;; Total initial energy
  (define e (bh-tree-total-energy eps energy-beta bs))
  
  ;; Timestep function.
  (define frac 0.5)
  (define timestep (max-force-timestep eps frac))
  
  ;; Single-step advancer.
  (define adv (lambda (bs dt) (gl3pt-advance eps advance-beta bs dt)))
  
  ;; Filter between steps---output roughly at output-interval.
  (define last-output (make-parameter 0))
  (define t (make-parameter 0))
  (define (filter bs dt)
    (t (+ (t) dt))
    (when (> (t) (+ (last-output) output-interval))
      (last-output (t))
      (let ((e-new (bh-tree-total-energy eps energy-beta bs)))
        (printf (string-append "At time ~a, energy is ~a, 10% lagrangian "
                               "radius is ~a, energy error is ~a~%")
                (t) e-new (lagrangian-radius bs 0.1) (abs (/ (- e e-new) 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 matrix-inverse)---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 self-contained). 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. 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 Street, Fifth Floor, Boston, MA 02110-1301 USA.
|#

(module base (lib "swindle.ss" "swindle")
  (provide (all-from-except (lib "swindle.ss" "swindle") + - * / add negate)
           (rename + mz:+) (rename - mz:-) (rename * mz:*) (rename / mz:/)
           (rename my-+ +) (rename my-- -) (rename my-* *) (rename my-/ /)
           (rename my-add add) sub mul div invert (rename my-negate negate))
  
  (define my-+
    (case-lambda
      (() 0)
      ((x) x)
      ((x y) (my-add x y))
      ((x y . ys)
       (apply my-+ (my-add x y) ys))))
  
  (define my--
    (case-lambda
      ((x) (my-negate x))
      ((x y) (sub x y))
      ((x . ys)
       (sub x (apply my-+ ys)))))
  
  (define my-*
    (case-lambda
      (() 1)
      ((x) x)
      ((x y) (mul x y))
      ((x y . ys)
       (apply my-* (mul x y) ys))))
  
  (define my-/
    (case-lambda
      ((x) (invert x))
      ((x y) (div x y))
      ((x . ys) (div x (apply my-* ys)))))
  
  (defgeneric (my-add x y))
  (defmethod (my-add (x ) (y ))
    (+ x y))
  
  (defgeneric (my-negate x))
  (defmethod (my-negate (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. 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 02110-1301, USA.

|#

(module matrix mzscheme
  (require (lib "42.ss" "srfi")
           (lib "contract.ss")
           (only "base.ss" struct-type->class add sub mul div invert negate defmethod
                  ))
  
  (define mz:+ +)
  (define mz:- -)
  (define mz:* *)
  (define mz:/ /)
  
  (define-struct matrix
    (rows cols elts) #f)
  
  (define  (struct-type->class struct:matrix))
  
  (provide struct:matrix :matrix matrix-ec *small* )
  
  (define (valid-matrix-row/c m) (and/c natural-number/c (between/c 0 (mz:- (matrix-rows m) 1))))
  (define (valid-matrix-col/c m) (and/c natural-number/c (between/c 0 (mz:- (matrix-cols m) 1))))
  
  (define (same-matrix-dimensions/c m)
    (flat-named-contract
     "same matrix dimensions"
     (let ((mc (matrix-cols m))
           (mr (matrix-rows m)))
       (lambda (m2)
         (and (= (matrix-cols m2) mc)
              (= (matrix-rows m2) mr))))))
  
  (define square-matrix/c
    (flat-named-contract
     "square matrix"
     (lambda (m) (= (matrix-rows m) (matrix-cols m)))))
  
  (provide/contract
   (rename my-make-matrix make-matrix (-> natural-number/c natural-number/c any/c matrix?))
   (rename my-matrix matrix (->r ((r natural-number/c)
                                  (c natural-number/c))
                                 elts (and/c list? (lambda (els) (= (length els) (mz:* r c))))
                                 matrix?))
   (matrix? (-> any/c boolean?))
   (matrix-ref (->r ((m matrix?)
                     (i (valid-matrix-row/c m))
                     (j (valid-matrix-col/c m)))
                    any))
   (matrix-set! (->r ((m matrix?)
                      (i (valid-matrix-row/c m))
                      (j (valid-matrix-col/c m))
                      (x any/c))
                     any))
   (matrix-mul (->r ((m1 matrix?)
                     (m2 (and/c matrix?
                                (flat-named-contract
                                 "compatible matrix dimensions"
                                 (lambda (m2) (= (matrix-cols m1) (matrix-rows m2)))))))
                    matrix?))
   (matrix-transpose (-> matrix? matrix?))
   (matrix-vector-mul (->r ((m matrix?)
                            (v (and/c vector?
                                      (flat-named-contract
                                       "compatible vector dimension"
                                       (lambda (v) (= (vector-length v) (matrix-cols m)))))))
                           vector?))
   (vector-matrix-mul (->r ((v vector?)
                            (m (and/c matrix?
                                      (flat-named-contract
                                       "compatible vector dimension"
                                       (lambda (m) (= (vector-length v) (matrix-rows m)))))))
                           vector?))
   (matrix-add (->r ((m1 matrix?)
                     (m2 (and/c matrix? (same-matrix-dimensions/c m1))))
                    matrix?))
   (matrix-sub (->r ((m1 matrix?)
                     (m2 (and/c matrix? (same-matrix-dimensions/c m1))))
                    matrix?))
   (matrix-scale (->r ((s (or/c number? matrix?))
                       (m (if (number? s) matrix? number?)))
                      matrix?))
   (matrix-norm (-> (and/c integer? (>/c 0)) (-> matrix? number?)))
   (matrix-infinity-norm (-> matrix? number?))
   (matrix-inverse (-> square-matrix/c square-matrix/c))
   (matrix-identity (-> natural-number/c matrix?))
   (matrix-copy (-> matrix? matrix?)))
  
  (define (my-make-matrix rows cols elt)
    (make-matrix rows cols (make-vector (mz:* rows cols) elt)))
  
  (define (matrix-ref m i j)
    (vector-ref (matrix-elts m) (mz:+ (mz:* i (matrix-cols m)) j)))
  
  (define (matrix-set! m i j x)
    (vector-set! (matrix-elts m) (mz:+ (mz:* i (matrix-cols m)) j) x))
  
  (define (my-matrix rows cols . elts)
    (make-matrix rows cols (list->vector elts)))
  
  (define-syntax matrix-ec
    (syntax-rules ()
      ((matrix-ec rrows ccols etc ...)
       (let ((rows rrows)
             (cols ccols))
         (make-matrix rows cols (vector-of-length-ec (mz:* rows cols) etc ...))))))
  
  (define-syntax :matrix
    (syntax-rules (index)
      ((:matrix cc var arg)
       (:vector cc var (matrix-elts arg)))
      ((:matrix cc var (index i j) arg)
       (:do cc
            (let ((m arg)
                  (rows #f)
                  (cols #f))
              (set! rows (matrix-rows m))
              (set! cols (matrix-cols 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 (matrix-mul m1 m2)
    (let ((rows (matrix-rows m1))
          (cols (matrix-cols m2))
          (n (matrix-cols m1)))
      (matrix-ec rows cols
        (:range i rows)
        (:range j cols)
        (sum-ec (:range k n)
          (* (matrix-ref m1 i k)
             (matrix-ref m2 k j))))))
  
  (define (matrix-transpose m)
    (matrix-ec (matrix-cols m) (matrix-rows m)
      (:range i (matrix-cols m))
      (:range j (matrix-rows m))
      (matrix-ref m j i)))
  
  (define (matrix-vector-mul m v)
    (vector-of-length-ec (matrix-rows m)
      (:range i (matrix-rows m))
      (sum-ec (:range j (matrix-cols m))
        (* (matrix-ref m i j) (vector-ref v j)))))
  
  (define (vector-matrix-mul v m)
    (vector-of-length-ec (matrix-cols m)
      (:range j (matrix-cols m))
      (sum-ec (:range i (matrix-rows m))
        (* (vector-ref v i) (matrix-ref m i j)))))
  
  (define (matrix-add m1 m2)
    (matrix-ec (matrix-rows m1) (matrix-cols m1)
      (:parallel (:matrix x1 m1)
                 (:matrix x2 m2))
      (+ x1 x2)))
  
  (define (matrix-sub m1 m2)
    (matrix-ec (matrix-rows m1) (matrix-cols m1)
      (:parallel (:matrix x1 m1)
                 (:matrix x2 m2))
      (- x1 x2)))
  
  (define (matrix-scale m s)
    (when (number? m)
      (let ((temp m))
        (set! m s)
        (set! s temp)))
    (matrix-ec (matrix-rows m) (matrix-cols m)
      (:matrix x m)
      (* x s)))
  
  (define ((matrix-norm n) m)
    (expt (sum-ec (:matrix x m)
            (expt (abs x) n))
          (/ n)))
  
  (define (matrix-infinity-norm m)
    (max-ec (:matrix x m) (abs x)))
  
  
  (define (scale-matrix-row! m i s)
    (do-ec (:range j (matrix-cols m))
      (matrix-set! m i j (* (matrix-ref m i j) s))))
  
  (define swap-matrix-rows!
    (case-lambda
      ((m i j pvec)
       (let ((oldi (vector-ref pvec i)))
         (vector-set! pvec i (vector-ref pvec j))
         (vector-set! pvec j oldi))
       (swap-matrix-rows! m i j))
      ((m i j)
       (do-ec (:range k (matrix-cols m))
         (let ((oldik (matrix-ref m i k)))
           (matrix-set! m i k (matrix-ref m j k))
           (matrix-set! m j k oldik))))))
  
  (define (add-scaled-matrix-row-to-row! m s i j)
    (do-ec (:range k (matrix-cols m))
      (matrix-set! m j k (+ (matrix-ref m j k) (* s (matrix-ref m i k))))))
  
  (define (pivot-row m j)
    (fold3-ec 'no-max (:range i j (matrix-rows m)) i (lambda (x) x) (lambda (i iold)
                                                                      (if (> (abs (matrix-ref m i j))
                                                                             (abs (matrix-ref m iold j)))
                                                                          i
                                                                        iold))))
  
  (define (matrix-pref m p i j)
    (matrix-ref m (vector-ref p i) j))
  
  (define (matrix-identity n)
    (matrix-ec n n (:range i n) (:range j n) (if (= i j) 1 0)))
  
  (define (matrix-copy m)
    (matrix-ec (matrix-rows m) (matrix-cols m) (:matrix x m) x))
  
  (define *small* (make-parameter 1.4901161193847656e-08 (lambda (sm) (if (not (>= sm 0))
                                                                          (error '*small* "Cannot set to negative number: ~a" sm)
                                                                          sm))))
  
  (define (matrix-inverse m-in)
    (let ((rows (matrix-rows m-in))
          (cols (matrix-cols m-in)))
      (let ((m (matrix-copy m-in))
            (id (matrix-identity rows))
            (pvec (vector-of-length-ec rows (:range i rows) i)))
        (do-ec (:range i rows)
          (begin
            (let ((pi (pivot-row m i)))
              (swap-matrix-rows! m i pi pvec)
              (swap-matrix-rows! id i pi))
            (let ((pelt (matrix-ref m i i)))
              (when (< (abs pelt) (*small*))
                (raise-mismatch-error 'matrix-inverse "Singular or nearly singular (by *small* parameter) matrix " (cons pelt m-in)))
              (scale-matrix-row! m i (/ pelt))
              (scale-matrix-row! id i (/ pelt))
              (do-ec (:range j rows) (not (= j i))
                (let ((s (- (matrix-ref m j i))))
                  (add-scaled-matrix-row-to-row! m s i j)
                  (add-scaled-matrix-row-to-row! id s i j))))))
        id)))
  
  (define (relative-size x norm)
    (abs (/ x norm)))
  
  (defmethod (add (m1 ) (m2 ))
    (matrix-add m1 m2))
  (defmethod (sub (m1 ) (m2 ))
    (matrix-sub m1 m2))
  (defmethod (mul (m ) (s ))
    (matrix-scale m s))
  (defmethod (mul (s ) (m ))
    (matrix-scale m s))
  (defmethod (mul (m ) (v ))
    (matrix-vector-mul m v))
  (defmethod (mul (v ) (m ))
    (vector-matrix-mul v m))
  (defmethod (mul (m1 ) (m2 ))
    (matrix-mul m1 m2))
  (defmethod (negate (m ))
    (matrix-scale m -1))
  (defmethod (invert (m ))
    (matrix-inverse m))
  (defmethod (div (v ) (m ))
    (matrix-vector-mul (matrix-inverse m) v))
  (defmethod (div (m ) (s ))
    (matrix-scale 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 :).

#| lorentz-transform.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 02110-1301, USA.

|#

(module lorentz-transform mzscheme
  (require "matrix.ss"
           (lib "42.ss" "srfi")
           (lib "contract.ss")
           (lib "math.ss")
           "simplex.ss"
           "simplex-map-utils.ss"
           (only "base.ss" struct-type->class defmethod  mul div invert)
           (all-except (planet "table.ss" ("soegaard" "galore.plt" 3))
                       map ormap andmap for-each values equal? union))
  
  (define-struct lorentz-transform-field
    (table) #f)
  
  (define-struct lorentz-transform-element
    (params inverse-params matrix inverse-matrix) #f)
  
  (define  (struct-type->class struct:lorentz-transform-field))
  (define  (struct-type->class struct:lorentz-transform-element))
  
  (provide eta lt-params/c struct:lorentz-transform-field 
           struct:lorentz-transform-element )
  
  (define lt-params/c (and/c vector?
                             (flat-named-contract
                              "6-element vector"
                              (lambda (v) (= (vector-length v) 6)))))
  
  (define zero-simplex/c (and/c simplex?
                                (flat-named-contract
                                 "simplex-dimension = 0"
                                 (lambda (s) (= (simplex-dimension s) 0)))))
  
  (define ltf-alist/c (listof (cons/c zero-simplex/c (or/c matrix? lt-params/c lorentz-transform-element?))))
  
  (provide/contract
   (matrix-eta-trace (-> matrix? number?))
   (lorentz-matrix (-> number? number? number?
                       number? number? number?
                       matrix?))
   (make-lorentz-matrix (-> lt-params/c
                            matrix?))
   (lorentz-matrix->params (-> matrix? lt-params/c))
   (Dlorentz-matrix (-> lt-params/c
                        (vectorof matrix?)))
   (lorentz-transform-field? (-> any/c boolean?))
   (lorentz-transform-element? (-> any/c boolean?))
   (params->lorentz-transform-element (-> lt-params/c lorentz-transform-element?))
   (matrix->lorentz-transform-element (-> matrix? lorentz-transform-element?))
   (lorentz-transform-element-params (-> lorentz-transform-element? lt-params/c))
   (lorentz-transform-element-inverse-params (-> lorentz-transform-element? lt-params/c))
   (lorentz-transform-element-matrix (-> lorentz-transform-element? matrix?))
   (lorentz-transform-element-inverse-matrix (-> lorentz-transform-element? matrix?))
   (lorentz-transform-element-compose (-> lorentz-transform-element?
                                          lorentz-transform-element?
                                          lorentz-transform-element?))
   (lorentz-transform-element-invert (-> lorentz-transform-element?
                                         lorentz-transform-element?))
   (lorentz-transform-element-transform-vector
    (-> lorentz-transform-element? vector? vector?))
   (rename my-make-lorentz-transform-field make-lorentz-transform-field
           (-> ltf-alist/c lorentz-transform-field?))
   (rename my-lorentz-transform-field lorentz-transform-field
           (->* () ltf-alist/c (lorentz-transform-field?)))
   (lorentz-transform-field-element (-> lorentz-transform-field? zero-simplex/c lorentz-transform-element?))
   (lorentz-transform-field-compose (-> lorentz-transform-field?
                                        lorentz-transform-field?
                                        lorentz-transform-field?))
   (lorentz-transform-field-fold (-> (-> zero-simplex/c lorentz-transform-element? any/c any)
                                     any/c
                                     lorentz-transform-field?
                                     any)))
  
  (define (matrix-eta-trace M)
    (+ (- (matrix-ref M 0 0))
       (matrix-ref M 1 1)
       (matrix-ref M 2 2)
       (matrix-ref M 3 3)))
  
  (define Rx-gen (matrix 4 4
                         0 0 0 0
                         0 0 0 0
                         0 0 0 -1
                         0 0 1 0))
  (define Ry-gen (matrix 4 4
                         0 0 0 0
                         0 0 0 -1
                         0 0 0 0
                         0 1 0 0))
  (define Rz-gen (matrix 4 4
                         0 0 0 0
                         0 0 -1 0
                         0 1 0 0
                         0 0 0 0))
  (define Bx-gen (matrix 4 4
                         0 -1 0 0
                         -1 0 0 0
                         0 0 0 0
                         0 0 0 0))
  (define By-gen (matrix 4 4
                         0 0 -1 0
                         0 0 0 0
                         -1 0 0 0
                         0 0 0 0))
  (define Bz-gen (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 Bx-gen By-gen Bz-gen Rx-gen Ry-gen Rz-gen))
  (define mat-funs (list Bx By Bz Rx Ry Rz))
  (define rgens (reverse gens))
  (define rmat-funs (reverse mat-funs))
  
  (define (make-lorentz-matrix param-v)
    (let ((! vector-ref)
          (* matrix-mul)
          (p param-v))
      (* (Bx (! p 0))
         (* (By (! p 1))
            (* (Bz (! p 2))
               (* (Rx (! p 3))
                  (* (Ry (! p 4))
                     (Rz (! p 5)))))))))
  
  (define (lorentz-matrix vx vy vz x y z)
    (make-lorentz-matrix (vector vx vy vz x y z)))
  
  (define (Dlorentz-matrix params)
    (let* ((lparams (vector->list params))
           (mats (map (lambda (f x) (f x)) mat-funs lparams))
           (rmats (reverse mats)))
      (vector-of-length-ec 6 (:range i 6)
        (fold-ec (matrix-identity 4) (:parallel (:list M rmats)
                                                (:list g rgens)
                                                (:range j 5 -1 -1))
          (if (= i j)
              (matrix-mul M g)
              M)
          matrix-mul))))
  
  (define (asinh z)
    (log (+ z
            (sqrt (+ (sqr z)
                     1)))))
  
  (define (acosh z)
    (log (+ z
            (* (sqrt (- z 1))
               (sqrt (+ z 1))))))
  
  (define (lorentz-matrix->params M)
    (let* ((sinh-vz (- (matrix-ref M 3 0)))
           (vz (asinh sinh-vz))
           (cosh-vz (cosh vz))
           (sinh-vy (- (/ (matrix-ref M 2 0)
                          cosh-vz)))
           (vy (asinh sinh-vy))
           (cosh-vy (cosh vy))
           (sinh-vx (- (/ (matrix-ref M 1 0)
                          (* cosh-vy cosh-vz))))
           (vx (asinh sinh-vx)))
      (let ((new-M (matrix-mul (Bz (- vz))
                               (matrix-mul (By (- vy))
                                           (matrix-mul (Bx (- vx)) M)))))
        (let* ((sin-y (- (matrix-ref new-M 1 3)))
               (y (asin sin-y))
               (cos-y (cos y))
               (sin-z (- (/ (matrix-ref new-M 1 2)
                            cos-y)))
               (z (asin sin-z))
               (cos-z (cos z))
               (cos-x (/ (matrix-ref new-M 3 3)
                         cos-y))
               (sin-x (/ (+ (matrix-ref new-M 3 2)
                            (* sin-y sin-z cos-x))
                         cos-z))
               (x (asin sin-x)))
          (vector vx vy vz x y z)))))
  
  (define (my-lorentz-transform-field . maps)
    (my-make-lorentz-transform-field maps))
  
  (define (params->lorentz-transform-element p)
    (let* ((M (make-lorentz-matrix p))
           (M-inv (matrix-inverse M))
           (p-inv (lorentz-matrix->params M-inv)))
      (make-lorentz-transform-element p p-inv M M-inv)))
  (define (matrix->lorentz-transform-element M)
    (let* ((M-inv (matrix-inverse M))
           (p (lorentz-matrix->params M))
           (p-inv (lorentz-matrix->params M-inv)))
      (make-lorentz-transform-element p p-inv M M-inv)))
  
  (define (my-make-lorentz-transform-field ltf-alist)
    (let ((elt-alist (map
                      (lambda (s-p-or-M)
                        (let ((p-or-M (cdr s-p-or-M)))
                          (cons (car s-p-or-M)
                                (cond
                                  ((vector? p-or-M)
                                   (params->lorentz-transform-element p-or-M))
                                  ((matrix? p-or-M)
                                   (matrix->lorentz-transform-element p-or-M))
                                  (else p-or-M)))))
                      ltf-alist)))
      (make-lorentz-transform-field (alist->ordered simplex-compare elt-alist))))
  
  (define (lorentz-transform-field-element ltf s)
    (lookup s (lorentz-transform-field-table ltf)
            (lambda ()
              (make-lorentz-transform-element
               (make-vector 6 0)
               (make-vector 6 0)
               (matrix-identity 4)
               (matrix-identity 4)))
            values))
  
  (define (lorentz-transform-element-compose elt1 elt2)
    (let* ((M1 (lorentz-transform-element-matrix elt1))
           (M2 (lorentz-transform-element-matrix elt2))
           (M (matrix-mul M1 M2))
           (p (lorentz-matrix->params M))
           (M-inv (matrix-inverse M))
           (p-inv (lorentz-matrix->params M-inv)))
      (make-lorentz-transform-element p p-inv M M-inv)))
  
  (defmethod (mul (e1 )
                  (e2 ))
    (lorentz-transform-element-compose e1 e2))
  
  (define (lorentz-transform-element-transform-vector
           elt v)
    (matrix-vector-mul (lorentz-transform-element-matrix elt) v))
  
  (defmethod (mul (e )
                  (v ))
    (lorentz-transform-element-transform-vector e v))
  
  (define (lorentz-transform-element-invert elt)
    (make-lorentz-transform-element
     (lorentz-transform-element-inverse-params elt)
     (lorentz-transform-element-params elt)
     (lorentz-transform-element-inverse-matrix elt)
     (lorentz-transform-element-matrix elt)))
  
  (defmethod (invert (e ))
    (lorentz-transform-element-invert e))
  (defmethod (div (v )
                  (e ))
    (matrix-vector-mul (lorentz-transform-element-inverse-matrix e) v))
  
  (define lorentz-transform-field-compose
    (let ((table-compose (simplex-map-lift-combiners
                          (lambda () (make-ordered simplex-compare))
                          insert
                          (lambda (sl el sr er accu)
                            (insert sl (lorentz-transform-element-compose el er) accu))
                          insert)))
      (lambda (ltf1 ltf2)
        (make-lorentz-transform-field
         (table-compose
          (lorentz-transform-field-table ltf1)
          (lorentz-transform-field-table ltf2))))))
  
  (defmethod (mul (ltf1 )
                  (ltf2 ))
    (lorentz-transform-field-compose ltf1 ltf2))
  
  (define (lorentz-transform-field-fold f start ltf)
    (fold f start (lorentz-transform-field-table 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 M-inv.
Compute the parameters associated with M-inv, p-inv.

This is a natural sequence of operations to perform on the group---it's just the group inverse function on the parameters. The timings for 10K repeats of this operation are:

Code Style	Total CPU Time	GC Time
Pure Scheme	51 sec	27 sec
Mixed Scheme/C	4 sec	2.8 sec

I think the large GC time for the Scheme/C code is probably due to the fact that each call goes through several contract checks, which probably dominate the run-time and allocation. Clearly, there's still an advantage to doing significant manipulations in C! I wish that it was practical to stay completely within PLT Scheme for doing this sort of work, but I'm afraid that it isn't. I suspect that Brad's implementation of choice (I assume), Gambit, would do much better on this sort of comparison.

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 idea---have 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.", linking with BLAS and LAPACK.

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 matrix-solve-many
    (get-ffi-obj 'dgesv_ *lapack*
                 (_fun (m b) ::
                       ((_ptr i _int) = (matrix-rows m))
                       ((_ptr i _int) = (matrix-cols b))
                       (_matrix = (matrix-copy m))
                       ((_ptr i _int) = (matrix-rows m))
                       (_u32vector o (matrix-rows m))
                       (x : _matrix = (matrix-copy b))
                       ((_ptr i _int) = (matrix-rows 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

Nothing to do with Scheme, but a really neat article here in the Washington Post. I'm really busy now, but I'm glad I took the time out to read this article. Found via the fascinating blog of Terrence Tao.

03 March 2007

Eager Comprehensions for Arrays

I just submitted a package to PlaneT for eager comprehensions (SRFI-42) for arrays (SRFI-25). It's a pretty short module, so I've also pasted it below.

#| array-ec.ss: SRFI-42 comprehensions for SRFI-25 arrays.
Copyright (C) 2007 Will M. Farr 

This 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  02110-1301  USA
|#

(module array-ec mzscheme
  (require (lib "42.ss" "srfi")
           (lib "25.ss" "srfi")
           (only (lib "1.ss" "srfi") fold))
  
  (provide (all-from (lib "42.ss" "srfi"))
           (all-from (lib "25.ss" "srfi"))
           array-ec :array)
  
  (define-syntax array-ec
    (syntax-rules ()
      ((array-ec shp expr expr2 ... final-expr)
       (let ((a (make-array shp)))
         (let ((a-rank (array-rank a)))
           (let ((sizes (list-ec (:range i a-rank) (- (array-end a i)
                                                      (array-start a i)))))
             (let ((index-map (lambda (i)
                                (apply values
                                       (car (fold
                                             (lambda (size indices-and-i)
                                               (let ((indices (car indices-and-i))
                                                     (i (cdr indices-and-i)))
                                                 (cons (cons (modulo i size) indices)
                                                       (quotient i size))))
                                             (cons '() i)
                                             (reverse sizes)))))))
               (let ((shared-a (share-array a (shape 0 (product-ec (:list size sizes) size)) index-map)))
                 (let ((i 0))
                   (do-ec expr expr2 ... (begin (array-set! shared-a i final-expr)
                                                (set! i (+ i 1))))
                   a)))))))))
  
  (define (make-index-generator arr)
    (let ((r (array-rank arr)))
      (let ((r-1 (- r 1))
            (lbs (vector-of-length-ec r (:range i r) (array-start arr i)))
            (ubs (vector-of-length-ec r (:range i r) (array-end arr i))))
        (let ((idx-v (vector-of-length-ec r (:vector lb lbs) lb)))
          (vector-set! idx-v r-1 (- (vector-ref idx-v r-1) 1))
          (lambda ()
            (vector-set! idx-v r-1 (+ (vector-ref idx-v r-1) 1))
            (let loop ((i r-1))
              (cond
                ((= i 0) (if (>= (vector-ref idx-v 0) (vector-ref ubs 0))
                             'done
                             idx-v))
                ((= (vector-ref idx-v i) (vector-ref ubs i))
                 (vector-set! idx-v i (vector-ref lbs i))
                 (let ((i-1 (- i 1)))
                   (vector-set! idx-v i-1 (+ (vector-ref idx-v i-1) 1))
                   (loop i-1)))
                (else idx-v))))))))
  
  (define-syntax :array
    (syntax-rules (index)
      ((:array cc x (index k0 ...) arr-expr)
       (:do cc
            (let ((arr arr-expr)
                  (gen #f))
              (set! gen (make-index-generator arr)))
            ((idx-v (gen)))
            (not (eq? idx-v 'done))
            (let ((i 0)
                  (x (array-ref arr idx-v))
                  (k0 #f)
                  ...)
              (begin (set! k0 (vector-ref idx-v i))
                     (set! i (+ i 1)))
              ...)
            #t
            ((gen))))
      ((:array cc x arr-expr)
       (:array cc x (index) arr-expr)))))

18 February 2007

Barnes-Hut Tree Code Submitted to PlaneT

I've just submitted some new code to PlaneT which implements Barnes-Hut Trees in PLT Scheme. It's surprisingly (to me) fast, creating a 100K-body tree in about 30 seconds on my old PowerBook G4 800 Mhz.

This is just another step in my code clean-up in preparation for comparing our algorithm to the standard leapfrog algorithm of cosmological codes.

N-Body Simulation Initial Conditions

I just packaged up some code for PlaneT which creates gravitational N-body initial conditions. I assume it'll appear in the next few days. It's not very extensive as of yet (only does cold and hot constant spatial density spherical distributions and the venerable Plummer model), but that's really all I've needed over the past year or so to test my algorithms.

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 Barnes-Hut-style cosmological simulations. The PlaneT package is a result of that consolidation.

My current plan is to gin up a prototype Barnes-Hut tree code in PLT Scheme (shouldn't take long---it'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 million-body 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.

Scheming