Figured it Out

The bit at the end of the last post where I was worried about the commutator (lie bracket in the lie algebra of vector fields, not vectors at the identity) of two left-invariant vector fields not being itself left-invariant was just a mistake. I shouldn't expect [extend(d/dtheta), extend(d/dphi)] = d/dpsi, but rather [extend(d/dtheta), extend(d/dphi)] = extend(d/dpsi)! If you do that, then it works out:

(test-case 
 "Extended tangent vectors are really left-invariant"
 (check-tuple-close? 
  1e-6
  ((vector-field->component-field 
    (- ((lie-algebra-bracket SO3) d/dphi d/dtheta)
       ((natural-extension SO3) d/dpsi))
    SO3-rectangular-chart)
   ((slot-ref SO3-rectangular-chart 'chiinv)
    (up 0.02345453 0.0349587 0.0435897)))
  (up 0 0 0)))

passes with flying colors. (In case the above code is not completely clear---irony of the century---it's computing the components of the vector field ([extend(d/dtheta),extend(d/dphi)]-extend(d/dpsi)) at some arbitrary point in SO3 and verifying that they're zero.

Lie Groups Mostly Working

It's working! (Mostly.) See:

;; Welcome to DrScheme, version 350.1-svn21jun2006.
;; Language: Swindle.

;; Require some modules 
(require "lie-group-SO3.ss")
(require "lie-groups.ss")
(require "manifolds.ss")
(require "tuples.ss")
;;

;;; Compute the structure constants for SO3 in the 
;;; x,y,z coordinate system.  Exactly what you 
;;; would expect.
(structure-constants SO3 SO3-rectangular-chart)
;#<down: elts=#(#<down: elts=#(#<up: elts=#(0 0 0)> 
;                              #<up: elts=#(0 0 1)> 
;                              #<up: elts=#(0 -1 0)>)> 
;               #<down: elts=#(#<up: elts=#(0 0 -1)> 
;                              #<up: elts=#(0 0 0)> 
;                              #<up: elts=#(1 0 0)>)> 
;               #<down: elts=#(#<up: elts=#(0 1 0)>
;                              #<up: elts=#(-1 0 0)> 
;                              #<up: elts=#(0 0 0)>)>)>
                                                                                                        
                                                                                                        
;; Name the euler angle coordinates
(define-named-coordinates (theta phi psi) 
  SO3-euler-angles-chart)

;; Note that coordinate vectors commute
((vector-field->component-field 
  (lie-bracket d/dtheta d/dphi) 
  SO3-euler-angles-chart) 
 (slot-ref SO3 'identity))
;; #<up: elts=#(0 0 0)>

;; While the extensions of coordinate vectors under 
;; left-multiplication do not.
((vector-field->component-field 
  ((lie-algebra-bracket SO3) d/dtheta d/dphi) 
  SO3-euler-angles-chart)
 (slot-ref SO3 'identity))
;; #<up: elts=#(0 0 1)>

[Edited: used to say that I didn't understand why [natural-extension(d/dtheta), natural-extension(d/dphi)] <> d/dpsi, but now I do.

Functional Differential Geometry in PLT Scheme

I've been working a bit lately with Gerry Sussman and Jack Wisdom to extend their software for functional differential geometry to handle Lie Groups (which, after all, are just manifolds). My thesis will probably have to do with treating General Relativity as an SO(3,1) gauge theory, so I'm particularly interested in the Special Orthogonal Lie Groups. (Jack, as a solar-system dynamicist, is also particularly interested in the Special Orthogonal groups because they represent rotations---in addition to being interested in GR as a gauge theory. Gerry is just interested in everything.)

Unfortunately, the software they have for doing this runs in MIT Scheme, and I own a PowerBook running OS X. It's not really a problem to SSH into a computer which runs their system, but it's nicer to have access to it on my own computer. So: I've coded up a bare-bones version of the scmutils/differential-geometry system for myself in PLT Scheme. It does no symbolic manipulation (only numerical calculations allowed), and doesn't have much of the nifty stuff that comes with scmutils proper, but it is able to compute on arbitrary manifolds and to handle Lie groups. I've been using SchemeUnit to run tests on the code as I go along, so there's a pretty good chance that it doesn't have major bugs.

I've posted my current darcs repository here; a darcs get http://web.mit.edu/farr/www/SchemeCode should get it for you, if you're interested. I haven't really written any documentation yet, but the test files should give some examples of how the functions are meant to be used. Comments are, of course, welcome at farr@mit.edu.

I'm not really sure why I'm posting this for the wider world right now (since it's so incomplete and "internal"). I hope someone finds it useful or interesting.

SRFI-78 for Bigloo

I've just posted an implementation of SRFI-78: Lightweight Testing for Bigloo. You can get it here. It leaves out the check-ec form because there's no implementation of SRFI-42 for Bigloo yet, but otherwise it's complete. The framework is very lightweight---it took about 10 minutes to port.

I want the framework so that I can test my new SRFI-4 implementation for Bigloo before I post it, so stay tuned for that! It feels good to be back in Scheme---much of my recent scientific coding has been in Ocaml (not that there's anything wrong with Ocaml, but I have missed Scheme).

Fun With Streams

Or: a bit of Haskell in PLT Scheme.

In a one-off version of the n-body integrator I've been working on (seemingly in perpetuity, but probably finishing soon with a couple of papers), I had occasion to iterate a function to convergence. (i.e. to compute (f (f (f ... (f x)))) until the result converges.) Well, it's pretty easy to write:

(define (iter-to-convergence f)
 (lambda (x)
   (let loop ((last (f x)))
     (let ((next (f last)))
       (if (= last next)
           next
           (loop next))))))

but where's the fun in that? It's much more fun to use streams for this.

So, using the streams module appended below (shamelessly stolen from SICP), I can write

(define (iter-to-convergence f)
 (lambda (x)
   (letrec ((vals
             (stream-cons (f x) (stream-map f vals))))
     (stream-take-converged vals))))

This idiom is very common in Haskell (and other lazy languages)---the canonical example is, of course, (define ones (stream-cons 1 ones)). The fact that this is so easy to add to scheme is a real testament to the language's power. What other language lets you (lazily) have your cake and (eventually) eat it, too?

Here's the code for the stream module (feel free to steal it if you want):

(module streams mzscheme
 (require (only (lib "1.ss" "srfi") any))

 (provide stream? stream-car stream-cdr stream-ref stream-cons
          (rename my-stream stream) stream-map stream-filter stream-take
          stream-converged? stream-take-converged)
  
 (define-struct stream (a b) #f)

 (define stream-car stream-a)

 (define (stream-cdr s)
   (force (stream-b s)))

 (define (stream-ref s n)
   (if (= n 0)
       (stream-car s)
       (stream-ref (stream-cdr s) (- n 1))))

 (define-syntax stream-cons
   (syntax-rules ()
     ((cons-stream car cdr)
      (make-stream car (delay cdr)))))

 (define-syntax my-stream
   (syntax-rules ()
     ((stream a)
      (stream-cons a ()))
     ((stream a b)
      (stream-cons a b))
     ((stream a b ...)
      (stream-cons a (my-stream b ...)))))      

 (define (stream-map f . ss)
   (stream-cons (apply f (map stream-car ss)) (apply stream-map f (map stream-cdr ss))))

 (define (stream-for-each f . ss)
   (apply f (map stream-car ss))
   (stream-for-each f (map stream-cdr ss)))

 (define (stream-filter test? s)
   (let ((s1 (stream-car s)))
     (if (test? s1)
         (stream-cons s1 (stream-filter test? (stream-cdr s)))
         (stream-filter test? (stream-cdr s)))))

 (define (stream-take n s)
   (if (= n 0)
       '()
       (cons (stream-car s) (stream-take (- n 1) (stream-cdr s)))))

 (define (stream-converged? s)
   (= (stream-car s) (stream-ref s 1)))

 (define (stream-take-converged s)
   (if (stream-converged? s)
       (stream-car s)
       (stream-take-converged (stream-cdr s)))))

Here's a more involved example of what you can do with this: the sieve of Eratosthenes (SICP 3.5.2).

(define (first-n-primes n)
 (stream-take
  n
  (letrec ((integers-from (lambda (n) (stream-cons n (integers-from (+ n 1)))))
           (primes (stream-cons 2 (stream-filter prime? (integers-from 3))))
           (prime? (lambda (x)
                     (let iter ((ps primes))
                       (let ((p (stream-car ps)))
                         (cond
                           ((= (remainder x p) 0) #f)
                           ((> (* p p) x) #t)
                           (else (iter (stream-cdr ps)))))))))
    primes)))