I think we’ve all been attracted to the Lorenz Attractor, but how many of us have had a closer look at an Ikeda map? Let me walk you through building an animated anti-aliased 2d renderer of mathematical functions (in this case Chaotic) in just 80 lines of Clojure! – Batteries included!
An Ikeda map is a discrete-time dynamic system. It looks similar to the well-known Lorenz attractor in that ‘vortexes’ appear and for certain values of u the behaviour is chaotic. The matematical formula is relatively simple:
One way to render the map is to select a number of random points in a coordinate system and feed these to the above equation as starting point. Then you iterate over them, producing new coordinates and this goes on forever. The only things that change from one rendering to the next is your starting points and the constant u. The position of the starting points doesn’t impact the rendering significantly, since they’re all sucked into the attractor points. For u=0.92 and after 1000 iterations, you get something similar to this:
Common for all these wonderful matemathical formulas, chaotic systems, fractals etc, is that we need some way to visualize them. There are 2 tools that I want right now,
- A 2d renderer
- A 3d renderer.
Let’s save #2 for a rainy day and get started on our 2D renderer.
Take a good hard look at the mathematical definition for an Ikeda map above, and then look at the Clojure translation:
(defn ikeda [x y u] (iterate (fn [[x_n y_n]] (let [t_n (- 0.4 (/ 6 (+ 1 (* x_n x_n) (* y_n y_n))))] [(inc (* u (- (* x_n (Math/cos t_n)) (* y_n (Math/sin t_n))))) (* u (+ (* x_n (Math/sin t_n)) (* y_n (Math/cos t_n ))))])) [x y]))
Start reading from the -let- statement. Firstly we define t_n and straight off you see the difference between infixed math (2 + 2) and prefixed math (+ 2 2). I remember a while ago I saw somebody trying to implement infixed math in Clojure, and I thought… Why? Once you’ve grown accustomed to the prefix notion, you don’t want to go back – You shouldn’t go back! From the -let- statement downwards, you’re seeing plain math. It’ll return a vector with the new coordinates ala [-0.1125 0.1123]. We pass that function to ‘iterate’ which conviently iterates it for us, seeding it with [x y]. In case iterate bothers you, compare with:
user> (iterate #(* 2 %) 1)) (1 2 4 8 16 32 64 128 256 512.......)
Now I’m basically done with the math. If I call (ikeda 10 10 0.902) I’ll get the ikeda map behaviour you saw in the picture above, for the point spawned at 10 10 – I’m now sitting with an infinite amount of iterations in my hand, just waiting to be computed. I think it’s worth noting, that once you get into this mind-set of building infinite sequences and pulling on them a world new world opens up in terms of programming, so it’s really worth embracing whole heartedly instead of trying to implement loops and what not.
Phew, that was the hard part. Moving on, we need to initialize our spawn points. One way to do this is simply to keep a collection of the lazy-streams:
(def spawns (for [i (range 200)] (let [x (nth (first axis-seqs) (rand-int (first dim-screen))) y (nth (last axis-seqs) (rand-int (last dim-screen)))] (atom (ikeda x y @u)))))
This means, that for every item in the collection ‘spawns’ I have an atom pointing to an infinite amount of iterations on that given point. Why an atom you ask? Well, because we’re holding 200 infinite sequences in our hand we will start to pile a lot of computations on the stack once we start realizing this data. To get around this, we need to give up the head of the sequence. To do that we render the first point and then set the atom to the ‘rest’ of it’s sequence. Basically like this:
user> [1 2 3 ... n] [1 2 3 ...] user> (rest [1 2 3 ... n]) (2 3 ... n) user> (rest (rest [1 2 3 ... n])) (3 ...)
Now the attentive reader will notice I’m calling the nth item out of axis-seqs. I thought it would be an easy read, if I generated a sequence which was exactly as wide and as high as my screen so that calling something like like nth 50 on the axis-seq, would tell me which coordinate the 50th pixel on my screen represented. The seqs are made like so:
(def axis-seqs [(iterate #(+ (/ (first dim-axis) (first dim-screen)) %) (first c-min)) (iterate #(+ (/ (last dim-axis) (last dim-screen)) %) (last c-min))])
In an imperative language you’d say ‘I know that each step along my math axis is width-axis / width-screen, so I’ll loop 800 times and increment my counter by that’, but in Clojure we define the mathematical formula for the increment, pass that to iterate and take out the number of items we want. Neat ? :)
With the math done and our spawn points generated, we move on to rendering. There’s only one problem with mathematical functions in terms of 2d rendering and that’s the axis. My screen goes from X: 0 -> 800, Y: 0->600, but my function is rendered in X: -20 -> 20, Y: -30 -> 30 for instance. So I need to map (-20,-30) -> (0, 0) and (20,30) to (800,600) as well as everything in between.
(defn >screencoordinate [coordinate] (map #(* (/ (- %1 %2) (- %3 %2)) %4) coordinate c-min c-max dim-screen))
Now I realize that this doesn’t read like plain english but if you understood the ikeda function, you ‘ll understand this. Calling (>screencodinate [5 -30]) will result in [444 0] if those were in fact my axis. The function ‘map’ takes a function as its first argument and a seq or series of seqs. On the first run it’ll pick out the first coordinate of each seq, lets say 5 and the first item of c-min which is -20, and so on. The math for the first step (X) then becomes :
So now that you can read Clojure like it was plain english it’s time to reap the benefits! How much would you have to change this code for work with 3 dimensions? or 4? or 5? You wouldn’t have to change one single character.
The UI stuff
Now, let’s code our UI system. This section demonstrates how to do some Java Interop and it’s more a demonstration of that, than anything else. Once it’s done however, it’s reusable.
First we need a function, which will render our data:
I’ll break my renderer into segments and commment on each directly below the code.
(defn draw-ikeda-map [#^Graphics canvas] (when (zero? @n) (doto canvas (.setColor Color/WHITE) (.fillRect 0 0 width-screen height-screen)))
We start out by checking ‘n’, which is our iteration counter. If it’s zero that means we’ve just begun our simulation and therefore need a clean background. I set the brush color to white and paint.In case you’re wondering about the #^Graphics in the argument list, that’s a type-hint. Because this function is performance critical I don’t want the compiler to do a reflection on this type, so I help it along.
Now I need to draw my spawn points:
(let [point-color (int (+ 155 (* (/ 100 iterations) (- @n iterations))))] (.setColor canvas (Color. point-color point-color point-color)) (doseq [spawnpoint spawns] (let [[x1 y1] (>screencoordinate (first @spawnpoint)) [x2 y2] (>screencoordinate (second @spawnpoint))] (.drawLine canvas x1 y1 x2 y2)))
I use a very customized home-made mathematical atrocity to determine which color to paint with. Basically the idea was to start out all black and slowly move towards white to distinguish the iterations. Once this is done I walk over my spawn-point collection, extracting the coordinates of the current iteration by calling (first ..) on our spawnpoint and then I get the next one by calling (second …). This is when the actual computation takes place. If the let statement looks crowded with brackets to you, that’s probably because you’ve haven’t been introduced to Clojures powerful destructuring semantics, but allow me to demonstrate:
user> (let [[a b c] [1 2 3]] b) 2 user> (let [[a b c & e] [1 2 3 4 5 6 7]] e) (4 5 6 7)
It’s an idiomatic way to split a collection into named variables without any language ceremony. In the first example I split [1 2 3] into [a b c] and return b, in the second I split 1 – 7 into [a b c] and what ever remains goes into e.
Finally, we need a HUD to see which iteration level we’re on and which value of ‘u’ we’re using:
(doto canvas (.setColor Color/WHITE) (.fillRect 0 3 100 30) (.setColor Color/BLACK) (.drawRect -5 2 101 31)) (.drawString canvas (format "u: %5f2" @u) 3 15) (.drawString canvas (format "iterations: %d" @n) 3 28)))
Again, it’s so straight forward it’s hardly worth explaining. I draw a white rectangle and an outline, and inside I print 2 strings. The doto macro expands each line into something like (.setColor canvas Color/WHITE).
Now there’s just one bit of logic left and that’s our animation thread. Basically we need a guy who keeps calling himself, asking spawns to drop the head, go to the next iteration, and ask for the canvas to be redrawn. I’ll do it step by step again:
(defn animate [surface] (doseq [point spawns] (swap! point next)) (swap! n inc)
First we walk over all points in our collection ‘spawn’ calling next on them. Since we’re dealing with atomic types we use ‘swap!’ to alter their value. Swap applies ‘next’ to the current value of the point and sets the atom to the resulting value. Remember the (rest [1 2 3]) example? That’s what’s going on.
Secondly our global atom ‘n’ which keeps track of our iteration level is inc(remented).
(when (>= @n iterations) (swap! u #(+ % (- 0.05 (float (/ (rand-int 100) 1000))))) (reset! n 0) (doseq [point spawns] (let [x (nth (first axis-seqs) (rand-int (first dim-screen))) y (nth (last axis-seqs) (rand-int (last dim-screen)))] (reset! point (ikeda x y @u)))))
Here we start by asking ‘Have we completed all the iterations we set out to do?’, if we have then we use a customized home-made mathematical atrocity to slightly modify our value of u, giving us a new simulation.
Then we call reset! on our global counter. Notice how swap! applies a function to the current value of the atom and sets the atom to the result, where reset! on the other hand, just sets a value.
With the globals reinitialized we need to reinitialize our spawn-points since they’re all computed from a certain value of u. This doseq loop is similar to the one at the top where we initially defined ‘spawns’.
(.repaint surface) (Thread/sleep 20) (when @running (recur surface)))
…We ask for the UI to be repainted while we artistically hold our breath. If our global ‘running’ is still true, we restart the thread by calling ‘recur’, and by implication if running is false, we stop.
The final interop!
Only one detail left. We need to fire up a JFrame, add a JPanel to it and override that Panels ‘render’ function so that we can supply our own. This is a practical detail, but I want to show it off anyway, because it’s saying a lot about Clojure’s quite impressive Interop capabilities:
(let [frame (JFrame. "Ikeda map") panel (doto (proxy [JPanel]  (paint [g] (render g))))] (doto frame (.add panel) .pack (.setSize (first dim-screen) (last dim-screen)) .show (.setDefaultCloseOperation JFrame/EXIT_ON_CLOSE)) (-> #(animate panel) Thread. .start))
Tada! With proxy I can easily override JPanels render method and use my own. With that out of the way, I init the components and launch the animator and that’s it. Not too bad considering you now have some very reuseable code for when you want to render your next complex function. The last line also shows you something which is quite unique to Clojure I believe. I call #(animate panel) to return animator-function (a ‘closure’). I pass this as the argument to construct a Thread and call start on the resulting instance – I doubt you’ll ever see such impressive Java Interop mixed with Closures in other languages!
To run this, type the following from a linux terminal:
$ java -cp clojure.jar clojure.main ikeda.clj
If you’re on Windows it’s basically the same approach, you just have to start by installing Linux, but afterwards you’ll also get this result in all its animated splendor:
Hope you enjoyed the read,
Source code for Ikeda.clj: here
Ps: Big thanks to the notorious CGrand who has been very helpful in many discussions, giving code-reviews, trimming lines etc etc.
It’s official: The commentators on this blog are fantastic!
First we had the situation in the Scala vs Clojure Concurrency post, where a commentator actually supplied the methodology which O’Reilley will use to correct the example in their book and now this…
#1: Zach Tellman has ported my source-code below to work with his own Clojure wrappers for the JOGL OpenGL library, so if you want to see this simulation mapped on a 3d rotatable disc, all you have to do is fetch his code: here, and follow the instructions on his Wiki to install the dependency.
Thanks guys for contributing!