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)))

Tidbits

My god! It's been about a month since I last posted (about PLT Scheme's new JIT). Just to prime the pump (hopefully more to follow more quickly):

• There's a very interesting discussion on the Ocaml mailing list about Software Transactional Memory for concurrent programming. In particular, see the referenced papers here and here. The basic idea is to wrap up all accesses to shared memory in a STM monad, which is eventually sequenced atomically. The neat thing about doing this is that the concurrent mechanisms compose; see the papers for how it all works out. Now, there's no reason another language couldn't do this, but Haskell is ideal because the pure semantics and the type system conspire to prevent any non-pure accesses from within the atomic section that are not to "registered" STM slots (i.e. the atomic sequencer won't thread parts of the IO monad through the computation---only STM monads are allowed). Pretty slick stuff.
• The nbody integrator paper I'm working on is coming along nicely (I'm a physicist, though you wouldn't know it to look at the blog); just some quick (i.e. a day or two long) simulations to run and then a lot of editing before it's ready for submission. Some time soon I may post a quick synopsis of the idea; for now you can read about what it's not (though similar) in Jack Wisdom's Symplectic Maps for the N-Body Problem (with Holman). If that tickles your fancy, have a look at Jack's Symplectic Correctors (also with Holman). If that doesn't make your jaw drop then you must not be much of a dynamicist!

PLT Gets JIT

According to this message to the PLT Scheme mailing list, PLT Scheme just got a JIT compiler. It only works on the non-3m (i.e. Boehm conservative collector) version of PLT Scheme, but it seems to provide some speedups. Subsequent messages suggest that they're working on a "proper" JIT compiler targeting LLVM which should provide even more speed. Very nice!

Update: I've been using the JIT version for the last day (my normal preference is for the 3m version), and it's really nice! I see speedup near a factor of 2 for much of my code (which is pretty loop-heavy, so maybe that helps get the extra speedups), and the whole experience of using DrScheme feels much more snappy!

Bitten By the Hygiene Bug

I've noticed that I haven't posted here for a while---been really busy working on the code for the paper on n-body integrators. I'm finally running some real "experiments" (the numerical kind), and am discovering the joy of tight binary orbits in the middle of a star cluster (hint: lots of small steps for the binary = lots of roundoff error over cluster timescales). Hopefully I've got that fixed up now, but you never know until it breaks :).

Anyway, here's a real life example of why hygiene matters. Consider the following code:

(define (vector-n-smallest vector n <)
   (let* ((store (vector-copy vector 0 (+ n 1))))
      (insert-sort! store <)
      (do-range (i (+ n 1) (vector-length vector))
         (vector-set! store n (vector-ref vector i))
         (insert-sort! store <))
      (vector-copy store 0 n)))

Hmm, my custom comparison function keeps complaining it's being fed integers rather than frobs! What's going on?

Turns out that I had defined do-range without macro hygiene, and it's expansion involves the symbol <. It's not rocket science, but it is a good example of why you should care about hygiene. (Note that this problem is "solved" in Common Lisp---I think---because the passed in comparison function lives in the value of the symbol <, while the CL::< numerical comparison lives in the symbol-function of the symbol <, so you have to use a funcall to get the custom comparison while (< ...) does the usual CL thing. Hygiene, on the other hand, just does the right thing.)

Bigloo 2.7a Released

In case you hadn't noticed: Bigloo 2.7a was released yesterday (or maybe over the weekend). You can get my OS X package for it if you don't like to compile things yourself.

Continuations in Unfriendly Environments

I'm sure many people already have seen the paper (and probably attended the conference), but have a look at this paper from Greg Pettyjohn, John Clements, Joe Marshall, Shriram Krishnamurthi, and Matthias Felleisen! It describes a technique that can be used to implement continuations in unfriendly environments where stack inspection is not available (i.e. portable ANSI C, JVM, .NET, etc.) The great thing about the technique (as opposed to, say, the CPS conversion in Chicken) is that most functions can still use the standard calling conventions and runtime stack of the host language or environment.

Here's my (probably weak) understanding of how it works: when a continuation is captured a chain of exceptions gets propagated up the stack. At each "stop" on the way up the stack, an exception handler cons-es up a closure which will continue the function whose frame is active at that point. When the exception reaches the top level of the stack, the frames are packaged up as the continuation. Activating the continuation simply calls the closures in order to resume the saved computation where it left off. All the implementation needs to provide is some way to unwind the stack, stopping at each frame along the way. In C, it would be a longjmp/setjmp; in Java it would be catch/throw; in .NET catch/throw.

There is an example with more code (and less math) here, too.

Exciting and Depressing

Ever have something happen that's terribly exciting, as long as you don't think about it too much? I just had a mental breakthrough with my integrators: I realized a way to think about what's going on that makes everything perfectly clear. Then I implemented this in about 174 lines of Bigloo code. "That's great!" you say. I agree completely. Unfortunately, if I can't quite remain lost in the moment, I think of two mildly depressing things:

1. Why didn't I think of this sooner? I've spent the last few months playing around with these things, and then I figure out that pretty much everything I was doing before was wrong in the space of two days---and it's (in retrospect) obvious!
2. Someone else has to have had this thought! I'm almost scared to do a literature search tomorrow (time for bed now)---I'm sure I'll find half a dozen papers describing this technique and explaining why it's not better than what everyone else has been doing, etc (worst case) or that other people have had this thought and it's got promise for future implementations (much better, but still a bit of a bummer).

However, I'm still terribly excited! 174 lines!! It's so simple!!! This is great!!!!

Shared State in PLT Scheme's Swindle

I just spent about six hours debugging a version of my integrator that runs in PLT Scheme (got tired of not having the class creation options available at runtime in Bigloo) because of a shared state problem. If you define a Swindle class with the slot option :initvalue <value> then <value> is only evaluated once; every object of this class shares the same value. Not good if you want to have mutable state! The thing to use is :initializer <thunk>. I wish there was an :initform slot option, like in the real CLOS.

SRFI-32 for Bigloo

I needed some fast vector sorting routines for my n-body integrator. To that end, I've packaged up the withdrawn SRFI 32 for Bigloo. The SRFI is technically withdrawn, but even half of a half-standard is better than none (and this way I don't have to delve into sorting myself---I'm sure Olin Shivers is way better at that than I am).

Wolfram Tones

Here's an idea sure to appeal to nerds everywhere (including me): cellphone ring tones generated by cellular automata. Though it's a bit hard to stomach how much Wolfram exaggerates his contributions to the science of complexity and the study of cellular automata and related computational algorithms, it's still an interesting idea. Discovered through this New York Times article.

In addition to the unpleasantness of enriching someone who seems---shall we say---reluctant to acknowledge other contributions to the study of cellular automata, I am too cheap to just drop $2 on a cellphone ring tone. However, I do think my custom tone is pretty neat.

Srfi-4 Vectors in Bigloo

I just thought I'd drop a quick note about srfi-4-like vectors in Bigloo. For those who don't know, srfi-4 specifies an interface for homogeneous numeric vector types. For those of us who do numerical computing, this is a big win---boxing a double, for example, every time you reference or store into a vector is just impossibly inefficient if you want to get anything numerical done. But, Bigloo doesn't have an implementation of srfi-4---so how is it that I'm such an enthusiastic user?

Here's how: Bigloo allows for the definition of typed vectors. Typed vectors are vectors that only hold one of the native Bigloo types. (Unfortunately, this cannot be done for user-created types like classes.) For example, the module clause

(module f64vector
   (type (tvector f64vector (double)))
   (eval (export-all)))

Creates homogeneous vectors of doubles labeled by the type ::f64vector and functions to access these vectors. Nearly the entire srfi-4 interface is created; the only warts are

1. The make-f64vector procedure requires two arguments; the initializing value is not optional.
2. The f64vector procedure is not created.

Note that the created ref and set! procedures are typed correctly, which helps immensely in the type analysis of the remainder of your program (I hardly ever end up putting in typing information except at procedure arguments).

In fact, the Bigloo source distribution contains an example of using tvectors, but it's not a widely known fact, so I thought I would mention it here for posterity's sake.

By the way, supposedly the type system in Bigloo can identify when vectors are used homogeneously, and automatically convert them to tvectors. That's nice in theory, but in practice it only works on vectors that are used in only a few procedures within the same module, and even then it has trouble if you use the vectors in some strange way or forget to explicitly initialize them in the creation statement. For things like numerical computing, where the f64vectors can get passed around through the entire program it's hopeless---just use a tvector.

Economic Integration

Check this out (from an article in the New York Times):

A school board election will take place in October. While the board has continued to endorse economic integration, some supporters worry that that could change one day.

"It's not easy and it can be very contentious in the community," said Walter C. Sherlin, who retired two years ago as an associate superintendent. "Is it worth doing? Look at 91 percent at or above grade level. Look at 139 schools, all of them successful. I think the answer is obvious."

This is great! A public initiative backed by research into student performance that seems to be wildly successful (40%->80% grade-level scoring among black students over the last decade, 79%->91% overall). And, unlike race-based busing, this doesn't suppose that your performance is the result of your race (which is an internal factor that you have no control over) but rather your economic situation (which is an external factor over which you have some control). I hope that the parents of the 2.5% of children who are involuntarily bussed don't become such a vocal minority that this stops happening.

Unemployment

Just read an interesting paragraph in an essay about the internet bubble:

The prospect of technological leverage will of course raise the specter of unemployment. I'm surprised people still worry about this. After centuries of supposedly job-killing innovations, the number of jobs is within ten percent of the number of people who want them. This can't be a coincidence. There must be some kind of balancing mechanism.

From Paul Graham. He's probably right, but it seems possible that the mechanism is people bitching about being unemployed and holding society effectively at ransom to find jobs for them to do. (i.e. a country that has 50% unemployment usually has social unrest, too, and some sort of reforms are proposed so that the rich people and companies can get back to making money without that pesky mob at the front gate.) Probably there's someone who studies this; it would be interesting to know whether my idea is right.

The Flu

Not that I have it (or usually get it, being young and in good shape). However, this article in the New York Times today points out that vaccines are "apparently ineffective" in preventing the flu for elderly people:

The fact that the vaccine study showed that inoculations have had only a modest effect in the elderly is particularly worrisome, because this is a group that tends to suffer high rates of complications and deaths from the disease and vaccination is the standard practice. In people over 65, the vaccines "are apparently ineffective" in the prevention of influenza, pneumonia and hospital admissions, although they did reduce deaths from pneumonia a bit, by "up to 30 percent," the study says.

"What you see is that marketing rules the response to influenza, and scientific evidence comes fourth or fifth," Dr. Jefferson said. "Vaccines may have a role, but they appear to have a modest effect. The best strategy to prevent the illness is to wash your hands."

Being young, I would like to see mention of the effectiveness of the vaccine in younger people (i.e. whether it's worth it for me to get one).

On a related note, how do people who don't believe in evolution at all (i.e. something made the world exactly as it is now) cope with the idea of viruses sharing genetic material and becoming resistant?

Bigloo Finalization Library

I've just posted a library which allows finalization of Bigloo objects before they're reclaimed by the GC. You can view the documentation and download the library here.

Manuel provided the hint to look at doc/README.finalizer, but I didn't like the way finalization was done in the example there, and it seemed a bit out of date. So, I wrapped up my own finalization library. I worry that I may have gone overboard (both in the documentation and in the care I've taken with license), but I have had a lot of trouble using poorly-documented libraries in the past, so hopefully this will be helpful to more people because I've spent some time writing out good documentation.

Bigloo 2.7a Beta Mac OS X Package

I've just posted a Mac OS X installer for the latest (Sept 20) beta of Bigloo. You can read more about it here. It's been quite a lot of work to compile Bigloo properly on OS X (which you can read about on the Bigloo mailing list), and Manuel has been a huge help. Many of these things are happening automatically now, but I figured that, since I've already put in the effort, I might as well make things as easy as possible for those following. Instead of putting a few #defines into some header files, typing ./configure --lots-of-options and then compiling, you can now just grab the disk image and Bigloo is only 4 clicks away. Enjoy!

Bigloo Internals, More Dynamic Linking and Shark...

without much Scheme-y stuff in sight.

Over the past few days Manuel Serrano has been helping me improve the speed of my n-body integrator by changing the behavior of Bigloo's threading when run in single-threaded mode on OS X. There is a function in runtime/Clib/cthread.c called denv_get which, I assume, gets some representation of the current dynamic environment. Since the cthread.c functions are really stubs for a more complete thread implementation denv_get is really short---it compiles to about 5 instructions (including prelude and post-lude) on my machine. Unfortunately, in the beta version of 2.7a from 18Aug denv_get was installed as the function to get the dynamic environment through a function pointer, preventing its inlining by the C compiler. My n-body code was spending upwards of 50% of its time in denv_get (remember: 5 instructions!). Manuel had me insert #define BGL_DYNAMIC_ENV_ALTERNATE 1 in runtime/Include/bigloo.h from the newest beta of 2.7a (12Sept)¹. This apparently allowed for inlining of denv_get in single-threaded mode---presto: faster by a factor of two.

I'm mentioning this for two reasons: 1. Manuel is really cool to take so much time to help me with my performance issue. 2. If you use OS X and Bigloo 2.7a beta, you should definitely get the 12Sept version and insert #define BGL_DYNAMIC_ENV_ALTERNATE 1 into runtime/Include/bigloo.h!

In other news, Maunel has used my patch fixing a problem with the build script which was preventing the building of the compiler with shared libraries even when the --sharedcompiler=yes option was provided to configure. This makes it much simpler to build a bigloo linked against dynamic libraries on OS X. All you have to do is make sure to provide --sharedcompiler=yes (and, if you like, --sharedbde=yes) to configure. Then edit the Makefile.config, changing the following variables:

• LD=gcc
• LDFLAGS=-dynamiclib -single_module
• LDLIBS=-lbigloogc-2.7a
• SHAREDSUFFIX=dylib

That's it---no restarting the make in the middle after changing options or anything. (I can't remember whether this works for building the Development Environment, but I do remember that if it fails it's pretty clear what to fix.)

Finally, I'd just like to put in a plug for SHARK, Apple's profiler that comes with OS X. You can profile a running process, obtaining a call graph with who-calls-who and elapsed time in each call path, see the memory usage pattern (and page faults), look at the disassembled object code with color-coded hot spots and pipeline stalls, and get helpful advice like using the reciprocal square root estimate assembly instruction rather than the math library sqrt when you don't need super accuracy. I can't imagine the amount of effort it would have been to do the profiling I've done today and yesterday with gprof! Bigloo works quite well with SHARK, since it outputs relatively idiomatic C code (though you might have to de-mangle function names with bigloo-demangle); other compilers like Gambit-C (C used as assembler) or Chicken (C used in CPS) wouldn't give you code that SHARK would help with much.

Eventually I'll get back to Scheme-y things (I haven't written a good macro in days!), but I thought it would be good to get down some of the things I've been doing lately.

1: You can get the 12Sept beta of 2.7a from INRIA's FTP Site or its mirror at U. Nice.

Hybrid Cars and Pork

In today's New York Times:

"I would like to get more hybrids out of our system because I do think it's something that is here to stay," William Clay Ford Jr., chairman and chief executive of Ford Motor, told reporters in Detroit on Aug. 23. Last fall, Ford introduced a hybrid version of its Escape sport utility vehicle, the first of several hybrids planned.

At the same time, the energy bill signed by President Bush on Aug. 8 effectively gave a break to American manufacturers by extending what could be a tax credit of as much as $3,400 per car to purchasers of the first 60,000 hybrids sold by a company. The credit phases out after that. Toyota sold more than 60,000 hybrids in the first six months of this year, so the tax law seems intended to help General Motors and Ford.

[...]

The new tax law "seems to benefit those who haven't done anything in the area of hybrids until now," Christopher Richter, a CLSA Asia-Pacific Markets auto analyst, said in Tokyo. "And it seems to penalize those who have been pioneers."

No kidding! And note that Ford has leased its hybrid technology from Toyota---no R&D costs for them. The only good news on the pork-ridden government front is that California is going to make a distinction between hybrids based on their environmental quality:

In early August, California regulators started to distinguish between fuel-efficient hybrids and "muscle" hybrids, the high-powered versions that save little gasoline. Of the seven hybrid models now on sale in the United States, owners of the Honda Civic, Honda Insight and Toyota Prius can qualify for decals allowing them to drive alone, rather than with two or more passengers, in highway commuter lanes. That reward is not extended to four hybrids not rated as exceptional energy savers: Honda Accord, Ford Escape, Toyota Highlander and Lexus RX 400h.