Mandelbrot in Scala

With my growing frustration with Matlab, I’ve been looking for a while for a language that was

  1. Garbage collected
  2. Good notation for numerical computation
  3. Fast enough for numerical computation

After a long search, I think I’ve finally found my home in Scala. Today I did a very trivial, self-contained computation of some images of the Mandelbrot set to test it out.

When learning Scala, I couldn’t find enough examples of numerical stuff, so I thought I would post this here for those potentially interested. Here’s a first attempt:

import java.io._
import java.lang.Math

object mandelbrot {
  // tiny complex number class including syntactic sugar for basic operations
  class Complex(val a: Double, val b: Double){
    // represents the complex number a + b*i
    def +(that: Complex) = new Complex(this.a+that.a,this.b+that.b)
    def *(that: Complex) = new Complex(this.a*that.a-this.b*that.b,this.a*that.b+that.a*this.b)
    def abs() = Math.sqrt(this.a*this.a + this.b*this.b)
  }

  def run(n: Int, level: Int) : Unit = {
    val out = new FileOutputStream("scalaimage.pgm")
    out.write(("P5\n"+n+" "+n+"\n255\n").getBytes())

    for {y0 <- 0 until n
	 x0 <- 0 until n }{
	   // y0 and x0 are in pixel integers
	   // x  and y  are real number coordinates
	   val x = -2.0 + x0*3.0/n
	   val y = -1.5 + y0*3.0/n

	   var z = new Complex(0,0)
	   var c = new Complex(x,y)
	   for(i <- 0 until level)
	     z = z*z + c 

	   if (z.abs < 2)
	     out.write(0);
	   else
	     out.write(255);
	 }
    out.close()
  }

  def main(args: Array[String]) {
    run(Integer.parseInt(args(1)), Integer.parseInt(args(0)))
  }
}

Not bad, right? Notice how painlessly we create the Complex class, along with the natural syntax for manipulating the numbers. What’s more, the result is very fast. For a 2048 by 2048 image, I get results in less than 30 seconds. (Just consider running this algorithm in Matlab. Ha!)

Save this to a file “mandelbrot1.scala”, and compile and run with:

scalac mandelbrot1.scala
scala mandelbrot 100 2048

The first argument is many times to try the mandelbrot iteration, and the second argument is how big of an image to output. The result is:

scalaimage11

(I manually converted the .pgm to .png to save space– click to get the full resolution version.)

A slightly more complex version follows. I added notation for the operators *= and +=, which speeds things up by about 10%. I also now color the pixels by how many iterations it takes the mandelbrot iterations to diverge.

import java.io._
import java.lang.Math

object mandelbrot {
  // tiny complex number class including syntactic sugar for basic operations
  class Complex(var a: Double, var b: Double){
    // represents the complex number a + b*i
    def +(that: Complex) = new Complex(this.a+that.a,this.b+that.b)
    def *(that: Complex) = new Complex(this.a*that.a-this.b*that.b,this.a*that.b+that.a*this.b)
    def abs() = Math.sqrt(this.a*this.a + this.b*this.b)
    def *=(that: Complex) ={
      val newa = this.a*that.a-this.b*that.b
      this.b = this.a*that.b+that.a*this.b
      this.a = newa
      this
    }
    def +=(that: Complex)={
      this.a += that.a
      this.b += that.b
      this
    }
  }

  def run(n: Int, level: Int) : Unit = {
    val out = new FileOutputStream("scalaimage.pgm")
    out.write(("P5\n"+n+" "+n+"\n255\n").getBytes())

    for {y0 <- 0 until n
	 x0 <- 0 until n }{
	   // y0 and x0 are in pixel integers
	   // x  and y  are real number coordinates
	   val x = -2.0 + x0*3.0/n
	   val y = -1.5 + y0*3.0/n

	   var z = new Complex(0,0)
	   var c = new Complex(x,y)
	   var i = 0
	   do {
	     z *= z
	     z += c
	     i += 1
	   } while( z.abs < 2 && i < level)

	   if (z.abs < 2)
	     out.write(0);
	   else
	     out.write( (i*255.0/level).toInt );
	 }
    out.close()
  }

  def main(args: Array[String]){
    run(Integer.parseInt(args(1)), Integer.parseInt(args(0)))
  }
}

Compile and run as before. The output is:

scalaimage2

I’m not sure how many people are using Scala for numerical applications, but it looks very good so far. There are a few minor tricks that turn out to be incredibly convenient. For example, if you define any function func for some class a, in a normal language you would call it by

a.func(b).

Scala defines this, but also gives you automatically the notation

a func b.

Adding in Scala’s non-mandatory type declarations (as above, notice that x0 and x aren’t declared to be Int or Double) tricks for implicit conversion between types, and it looks like it is possible to recreate almost all of Matlab’s synactic sugar for matrix manipulation, while actually allowing reasonable speed for non-vectorizable code, and reasonable facilities for abstraction. I hope to soon be an ex-Matlab user.

The scalala library looks very promising, though I haven’t taken the time to install it yet, and it looks like absolutely no one is using it. (dramage, give us some documentation!)

Advertisements

Unbiased coinflips from biased coinflips

A old problem, due to von Neumann goes as follows:

You have a biased coin that produces heads (H) with probability p, and tails (T) with probability (1-p).You don’t know p. How can you use this coin to simulate an unbiased coin?

The next paragraph contains a solution, so if you want to solve the problem yourself, stop reading now!

von Neumann’s solution was as follows*:

  1. Flip the coin twice
  2. If the outcome is HT, output H
  3. If the outcome is TH, output T
  4. Otherwise, go to 1.

A nice solution, but you can see that you might need to flip the coin many times. In particular, the probability of getting either HH or TT in any particular round is p^2 + (1-p)^2, which could be really big for highly biased coins.

Is there a more efficient method?

Today I found a beautiful paper examining this question. The insight is that the von Neumann scheme is based on symmetry– picking pairs of output strings that have equal probability, then outputting heads for one and tails for the other.

You can draw an (infinitely large) tree, with branches corresponding to random coin flip outcomes, and leaf nodes for an output. The expected number of coin flips for an output is the expected depth in the tree that one reaches before outputting a result and quitting. Viewed in this light, one can think of ways to create more efficient algorithms. The paper contains trees illustrating this wonderfully.

* Incidentally, when starting a paragraph with the name of a person whose first name is not capitalized, should the first character of the paragraph be capitalized? Which convention takes precedence?