15 November 2006

The Paper is Out!

It's out! The paper of the century---get your copy before they're sold out at the newsstand. (If you're into N-body simulation methods, that is.)

Seriously: I'm very glad to have it out the door. The ball is now in ApJ's court, and I can get back to focusing on my thesis (for which I just submitted a title: "Numerical Relativity from a Gauge Theory Perspective").

Section 6.2.1 uses the OCaml code I posted about earlier for automatic differentiation to show that the integrator algorithm is as symplectic as we claim. This is something I'd like to see more physicists take up---particularly in the cosmological N-body simulations, where preserving the "coldness" of the dark matter phase-space distribution is very important. Everyone uses the leapfrog algorithm (for now---it might turn out that our algorithm is better) because it's symplectic at constant timesteps (and, though we're the first to realize it, with block-power-of-two timesteps---see Section 3, and numerical experiments in Sections 6.1 and 6.2.1), but they never measure whether this is really the case with their algorithms. If they were in the habit of doing this, I'm sure someone before us would have realized that block-power-of-two steps preserve symplecticity while ordinary adaptive steps break it (rather than assuming that any non-constant timestep breaks it).

Anyway, I'm excited to see what the reception of the physics community will be. I've got my fingers crossed....

07 November 2006

SRFI-42 Comprehensions for SRFI-4 Vectors

I spend quite a bit of time working with SRFI-4 homogeneous numeric vectors. Below is some code which adds some SRFI-4-vector generators and comprehensions to SRFI-42's Eager Comprehensions. I am impressed at the modularity in SRFI-42!

#|
Derived (almost verbatim) from the code for vector comprehensions in the SRFI-42 reference implementation.  The copyright on that code is:

Copyright (C) Sebastian Egner (2003). All Rights Reserved.

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Modifications are copyright (C) Will M. Farr (2006), goverened by the same conditions.  Feel free to email comments to .  
|#

(module srfi-4-comprehensions mzscheme
  (require (lib "42.ss" "srfi")
           (lib "4.ss" "srfi"))
  
  (provide :s8vector :u8vector :s16vector :u16vector :s32vector 
           :u32vector :s64vector :u64vector :f64vector :f32vector
           
           s8vector-ec u8vector-ec s16vector-ec u16vector-ec s32vector-ec 
           u32vector-ec s64vector-ec u64vector-ec f64vector-ec f32vector-ec
           
           s8vector-of-length-ec u8vector-of-length-ec s16vector-of-length-ec 
           u16vector-of-length-ec s32vector-of-length-ec 
           u32vector-of-length-ec s64vector-of-length-ec u64vector-of-length-ec 
           f64vector-of-length-ec f32vector-of-length-ec)
  
  (define-for-syntax symbol-append
    (case-lambda
      ((s) s)
      ((s1 s2 . ss)
       (apply symbol-append (string->symbol (string-append (symbol->string s1) (symbol->string s2))) ss))))
  
  (define-syntax make/prefix
    (lambda (stx)
      (syntax-case stx ()
        ((make-prefix-generator prefix)
         (let ((pre-sym (syntax-object->datum (syntax prefix))))
           (with-syntax ((vlength (datum->syntax-object (syntax prefix) (symbol-append pre-sym 'vector-length)))
                         (vref (datum->syntax-object (syntax prefix) (symbol-append pre-sym 'vector-ref)))
                         (vgen (datum->syntax-object (syntax prefix) (symbol-append ': pre-sym 'vector)))
                         (vfilter (datum->syntax-object (syntax prefix) (symbol-append 'ec-: pre-sym 'vector-filter)))
                         (vmake (datum->syntax-object (syntax prefix) (symbol-append pre-sym 'vector-make)))
                         (vset! (datum->syntax-object (syntax prefix) (symbol-append pre-sym 'vector-set!)))
                         (v->list (datum->syntax-object (syntax prefix) (symbol-append pre-sym 'vector->list)))
                         (list->v (datum->syntax-object (syntax prefix) (symbol-append 'list-> pre-sym 'vector)))
                         (v-ec (datum->syntax-object (syntax prefix) (symbol-append pre-sym 'vector-ec)))
                         (v-of-length-ec (datum->syntax-object (syntax prefix) (symbol-append pre-sym 'vector-of-length-ec))))
             (syntax 
              (begin  
                (define-syntax vgen
                  (syntax-rules (index)
                    ((vgen cc var arg)
                     (vgen cc var (index i) arg) )
                    ((vgen cc var (index i) arg)
                     (:do cc
                          (let ((vec arg) (len 0)) 
                            (set! len (vlength vec)))
                          ((i 0))
                          (< i len)
                          (let ((var (vref vec i))))
                          #t
                          ((+ i 1)) ))
                    ((vgen cc var (index i) arg1 arg2 arg (... ...))
                     (:parallel cc (vgen cc var arg1 arg2 arg (... ...)) (:integers i)) )
                    ((vgen cc var arg1 arg2 arg (... ...))
                     (:do cc
                          (let ((vec #f)
                                (len 0)
                                (vecs (vfilter (list arg1 arg2 arg (... ...)))) ))
                          ((k 0))
                          (if (< k len)
                              #t
                              (if (null? vecs)
                                  #f
                                  (begin (set! vec (car vecs))
                                         (set! vecs (cdr vecs))
                                         (set! len (vlength vec))
                                         (set! k 0)
                                         #t )))
                          (let ((var (vref vec k))))
                          #t
                          ((+ k 1)) ))))
                (define (vfilter vecs)
                  (if (null? vecs)
                      '()
                      (if (zero? (vlength (car vecs)))
                          (vfilter (cdr vecs))
                          (cons (car vecs) (vfilter (cdr vecs))) )))
                (define-syntax v-ec
                  (syntax-rules ()
                    ((v-ec etc1 etc (... ...))
                     (list->v (list-ec etc1 etc (... ...))) )))
                (define-syntax v-of-length-ec
                  (syntax-rules (nested)
                    ((v-of-length-ec k (nested q1 (... ...)) q etc1 etc (... ...))
                     (v-of-length-ec k (nested q1 (... ...) q) etc1 etc (... ...)) )
                    ((v-of-length-ec k q1 q2             etc1 etc (... ...))
                     (v-of-length-ec k (nested q1 q2)    etc1 etc (... ...)) )
                    ((v-of-length-ec k expression)
                     (v-of-length-ec k (nested) expression) )
                    ((v-of-length-ec k qualifier expression)
                     (let ((len k))
                       (let ((vec (vmake len))
                             (i 0) )
                         (do-ec qualifier
                                (if (< i len)
                                    (begin (vset! vec i expression)
                                           (set! i (+ i 1)) )
                                    (error "vector is too short for the comprehension") ))
                         (if (= i len)
                             vec
                             (error "vector is too long for the comprehension") ))))))))))))))

  (make/prefix s8)
  (make/prefix u8)
  (make/prefix s16)
  (make/prefix u16)
  (make/prefix s32)
  (make/prefix u32)
  (make/prefix s64)
  (make/prefix u64)
  (make/prefix f32)
  (make/prefix f64))