(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.
27 June 2006
26 June 2006
;; 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.
24 June 2006
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 email@example.com.
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.
20 June 2006
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).
08 June 2006
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)))