Sounds like a neat program. I'm on a laptop right now but i'll check it out later. The reason I am mailling is because you can use mencoder to convert a stream of image files into a video file. Indeed, it is pretty neat. I'd post an image, but I'm not sure whether

Ryan Dickie wrote: the other people on this list would appriciate a binary attachment. I'm hoping to make a DVD of various simulations - but that's kind of difficult when rendering full-size animations takes many hours! >_< Hence the request for optimisation help... ;-) Mencoder works on Linux, IrfanView + VirtualDub does it nicely on Windoze, I'm sure MacOS has something that can stitch PPM images together too. Use whatever you have on your platform. :-D

Andrew Coppin shares with us his Chaos program. I confirm that on my HP laptop it was faster than 15 minutes, but I won't speculate how to optimize it. I appreciated the elegance of overloading, the usage of Num classes, etc, which makes it more readable, although somewhat slower. Actually it took slightly less than 10 minutes for 300x300. I wonder whether making some constants global (such as dt) would change anything. For those who don't have patience to execute the program, I converted the ppms to a XVID coded AVI file. Thanks, Andrew http://users.info.unicaen.fr/~karczma/Work/Chaos0.avi What I didn't appreciate was the use of simple extrapolating Euler's method which for oscillating systems is known to be unstable, so the results of the simulation may be far from the reality. Well, one chaos is worth another one, and the sin is not as mortal as in the case of truly periodic systems, but it may be the cause that it is difficult to see the classical fractal structure of the attraction domains on the generated images. Try to use leapfrog Verlet... It will require to keep not only p and v for each frame, but also the acceleration, but the corrections are minor. Of course it will be slower, but then, why not increase dt? Jerzy Karczmarczuk

Andrew Coppin wrote: [snip]

Fact #1: I don't *know* of any other numerical integration algorithm. (I've heard of RK4, but it's too complicated for me to understand.)

higher-order runge-kutta algorithms typically have high integration accuracy, but may require multiple energy/force evaluations per integration step. also, they are generally not symplectic [1, 2], which is believed to be rather undesirable if you're looking for long-term stability in your simulation. stoermer/verlet/leap-frog is symplectic, and requires only one force evaluation per time-step. it is therefore _the_ default choice for numerical integration of the equations of motion. it is really quite simple [3, 4]. if you're _really_ looking for speed, you'll want to look at r-RESPA, which is a generalization of leapfrog that allows for multiple-timestep integration [5]. if the force evaluation is indeed the limiting step in your computation (which i'm not sure is the case), you'll certainly want to look into that. it basically boils down to splitting your interaction function into a (rugged) short-range part (which is zero beyond some distance) and which is evaluated every timestep, and a (smooth) long-range part, which is evaluated only every so many steps.

Fact #2: I have tried running the simulation with several different, non-comensurate time step values, and it always seems to produce the same output, so I'm reasonably confident there are no integration errors.

getting a good feeling for the maximum allowed timestep size is usually the first step to getting performance out of your program. as a rule of thumb, you'll want more than 10 steps per "period" of your fastest degree of freedom. you can get a feeling for that by finding the maximum second derivative of the potential energy function (which is going to be close to one of your magnets, i suppose. speaking of which -- does the energy diverge near the magnets? i think it does.). to validate the above, you'll want to look at energy conservation -- compute the total (potential and kinetic) energy of your system at the beginning and the end (for a given simulation length in terms of time, not in terms of timesteps) of your simulation, and compare to the starting energy. you'll find that there's a critical step size below which your energy stays constant, and above which it starts to diverge rather quickly. you might want to take that step size (and divide it by half if you're cautious). kind regards, v. [1] http://srv.chim.unifi.it/orac/MAN/node5.html [2] http://mitpress.mit.edu/SICM/book-Z-H-57.html#%_chap_5 [3] http://einstein.drexel.edu/courses/CompPhys/Integrators/leapfrog/ [4] http://shootout.alioth.debian.org/gp4/benchmark.php?test=nbody&lang=all [5] http://cat.inist.fr/?aModele=afficheN&cpsidt=15255925

Simon Peyton-Jones wrote:

| The source code has explicit monomorphic types all over it; I would | expect GHC to be able to optimise out any method calls. (OTOH, I'm not a | GHC expert... Simon? Don?)

Full monomorphisation is a whole-program optimisation and GHC isn't (yet) a whole-program compiler. Overloaded functions are indeed specialised for the types at which they are called *in the module they are defined*.

Simon

What I'm trying to say is, I would have expected, say, resultant :: [Vector2] -> Vector2 resultant = sum to run faster than resultant = sum In the latter case, we have resultant :: (Num n) => [n] -> n. I'm no expert, but I would expect that this requires a method lookup on each call, whereas with the more restrictive type I would expect the compiler to be able to just link directly to the correct method implementation. Am I mistaken?

| In the latter case, we have resultant :: (Num n) => [n] -> n. I'm no | expert, but I would expect that this requires a method lookup on each | call, whereas with the more restrictive type I would expect the compiler | to be able to just link directly to the correct method implementation. | Am I mistaken?

Probably not. Remember that sum is *itself* overloaded, so if resultant calls sum there is not much it can do except pass on the appropriate dictionary to sum.

Better would be to inline (and specialise) sum at this call site. But at the moment GHC doesn't specialise functions defined in *other* modules. That's another thing that could be improved (at the cost of risking constructing multiple copies of the same specialisation of 'sum'.

Right. I see what you're saying. Presumably changing it to resultant = foldl' (+) wouldn't help much either? Seems to me there's always a tradeoff to be made between CPU time, RAM usage, and code size - if not other factors too. :-S But Haskell seems to be fairly unique in that (it looks like) it's possible to transform the code quite aggressively without too much difficulty. You wouldn't ever want to write source code for six different versions of sum, but if a compiler can do it automatically that's another matter. ;-) I guess I was assuming that a function like sum is "simple enough" to get inlined at the call site - and hence possibly optimised at that site. Apparently not. By the way, is there some way of discovering what GHC has "really" done with your code, rather than what you "think" it's done? (Short of getting out a dissassembler that it...)

Simon Peyton-Jones wrote:

| I guess I was assuming that a function like sum is "simple enough" to | get inlined at the call site - and hence possibly optimised at that | site. Apparently not.

There is no reason in principle why not. It's just that GHC doesn't do it at the moment.

On one hand, it looks like something that will often be a win. On the other hand, how do you prevent a combinatorial explosion of specialised functions that hardly ever get called? Hmm... maybe this is why I'm not a compiler designer? ;-)

| By the way, is there some way of discovering what | GHC has "really" done with your code, rather than what you "think" it's | done?

Yes: -ddump-simpl.

Quoting the GHC manual: "HACKER TERRITORY. HACKER TERRITORY. (You were warned.)" Ah, that made me smile... Boy, does this thing generate a lot of output! Ah well, I was warned. ;-) I suppose to really make any sense of this lot, I would need to have a very deep understanding of the low-level internal workings of GHC. (I have read a few papers about GHC, but most of them were far too technical for me to comprehend.) Well anyway, I see why adding -O2 makes it run so much faster. When I add that flag, the output goes from (+) :: Vector3 -> Vector3 -> Vector3 to something like (+) :: GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> (# GHC.Float.Double, GHC.Float.Double, GHC.Float.Double #) and the body went from being an incomprehensible tangle of deeply nested case expressions to being a recognisable expression involving GHC.Prim.+##. (Anybody know what the difference between GHC.Prim.Double# and GHC.Float.Double is?) This in spite of the fact that the definition is actually (+) = vzip (+) vzip f (Vector3 x0 y0 z0) (Vector3 x1 y1 z1) = Vector2 (f x0 x1) (f y0 y1) (f z0 z1) So GHC has "unrolled" several levels of indirection here, and replaced an algebraic datatype with a bunch of primitive arguments. In other words, GHC has done what I wanted it to do. :-D There is an absolutely *huge* chunk of code for the derived instance of Read Colour. (Why did I add that again? Hmm, I think I just included the standard library version of my code. Ah well, it's not hurting anything.) Also moderately puzzled by things like "Str: DmdType U(U(L)U(L)U(L))U(U(L)U(L)U(L))m". Surely that's hacker territory if ever I saw it. ;-) Ah well, I doubt I'm going to come up with any new ideas for how to make my code go faster, but it's mildly entertaining wading through over 200 KB of textual output trying to guess what it means.

| Ah well, I doubt I'm going to come up with any new ideas for how to make | my code go faster, but it's mildly entertaining wading through over 200 | KB of textual output trying to guess what it means. Ha ha. You did say that you wanted to know what GHC *really* does :-) Seriously, there is mileage in someone working on the core-language pretty-printer to make it less intimidating. It was never really designed for wide consumption, and could be a lot better. Suppressing more qualified names, for example. Simon

| I wonder... is there a high-level document anywhere which explains what | GHC actually does with your code when you compile it? A lot, actually: http://hackage.haskell.org/trac/ghc/wiki/Commentary

On 07/05/07, Andrew Coppin wrote:

(Anybody know what the difference between GHC.Prim.Double# and GHC.Float.Double is?)

It's the difference between unboxed and boxed types. A boxed type's representation is in fact a pointer to the unboxed type (I think), so that a Double would be internally represented to a pointer to a Double#. As a consequence, a Double can be _|_, because this can be represented by a null pointer. So such luck with unboxed types. So working with unboxed types is quicker and consumes less memory, but don't use them in any kind of high level because the lack of a _|_ will bite you sooner or later. -- -David House, dmhouse@gmail.com

On 07/05/07, David House wrote:

represented by a null pointer. So such luck with unboxed types. So working with unboxed types is quicker and consumes less memory, but don't use them in any kind of high level because the lack of a _|_ will bite you sooner or later.

I hesitate to say that... it might come back and bite you in the _|_? I'll get me coat... D.

Try adding strictness annotations to all the components of all your data structures (i.e. put a ! before the type). Not all of the need it, but I doubt any need to be lazy either. Probably the reason quant8 seems to be taking so much time is that it is where a lot of stuff finally gets forced. Certainly, for things that are "primitive" like Colour and Vector you want the components to be strict, in general. Yes, originally the profile was showing quant8 taking something absurd

like 80% of the CPU time. When I changed the framebuffer to an IOUArray, the time spent in quant8 dropped *drastically*. (Because now the framebuffer is strict, and that's forcing the evaluation sooner.) I could certainly try making vectors, colours and arrays strict and see if that does something... (Thinking about it, the colour computation has a square root in it, and I bet that doesn't get forced until it hits quant8... Square root is an expensive operation on currentl hardware isn't it?)

Also, you have many many superfluous parentheses and use a different naming convention from representative Haskell code (namely camelCase). This is a pet hate of mine. NamesLikeThis are fine. names_like_this are fine too. But for the love of God, namesLikeThis just looks stupid and annoying! So I generally use camel case for stuff which has to start uppercase, and underscores for stuff that has to start lowercase. It's a system, and it works. Unfortunately it's not the standard convention in Haskell. (And I doubt I will convince anybody to change it...)