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

This entry was posted in Uncategorized and tagged , , . Bookmark the permalink.

28 Responses to Mandelbrot in Scala

  1. pizer says:

    Yes, Matlab is slow when you need to touch every scalar by yourself (as opposed to utilizing its vector and matrix operations). Apart from that I really like it for prototyping & testing algorithms.

    The number of people using Scala for numerical applications is probably very close to zero. Of course, I understand that every language has its use (somewhere). But since you mentioned number chrunching as application I HAVE to suggest — forgive me for sounding like a fan-boy — C++ as solution. You don’t need garbage collection for number chrunching at all. Seriously. With C++ you get operator overloading, can use abstractions without memory/speed penalties, and have a convenient language feature that helps you a great deal with managing resources like memory: RAII.

    Cheers!
    P

  2. James Iry says:

    With a little implicit magic, you can even make complex numbers look nice at a cost of extra object creation.

    object Complex {
    val i = new Complex(0,1)
    implicit def double2Complex(real : double) = new Complex(real, 0)
    }

    val c1 = 4.0 + 5.0 * i
    val c2 = -1.0 + 2.7 * i

    There’s been some talk from Martin Odersky about optimizing conversions like that so that no actual object creation takes place for conversions like that when it’s not necessary.

  3. justindomke says:

    @pizer- I should have probably added a 4th criteria: “not C++”. I do really like Matlab for the things it is natural for, but lately I’ve been slamming against its limitations. Because of this, I’ve been working in C++ and I… I really want garbage collection. Maybe that’s because I am uninformed. For example, I wrote my own buggy little matrix library with matlabish notation. I can’t possibly see how to do that with out using ‘new’. (Although as a consumer of the library I rarely need to manage memory myself, which might be your point.) I also ended up having to conservatively make copies of stuff, whereas with garbage collection I would never have to do this.

    @james- thanks!

  4. pizer says:

    I believe Matlab uses ‘copy on write’. If you want that you’ll need to invoke ‘new’ somewhere and use reference counting a la std::tr1::shared_ptr. In any case std::vector could be used to manage the life-time of the matrix’ coefficients. Since std::vector can be a class member you’re all good.

    Anyway, it was merely a suggestion. I guess it depends on your definition of “fast enough”. I have to handle big nonlinear optimization problems. It can never be fast enough. Also, memory consumption is an issue. Btw: I use Boost.uBLAS for linear algebra.

    Cheers!
    P.

  5. Break iterration c.abs > 2 because then it diverge. This will turbocharge the implemenation, won’t it?

    for(i 2)

    And then why don’t get rid of the square root at all?

    def absSqr() = this.a*this.a + this.b*this.b

    and test for c.absSqr < 4 instead of c.abs < 2

  6. justindomke says:

    @Philipp

    Good point about changing abs() to absSqr()! (I should have thought of that.) I tried it, but I got a pretty small speed improvement of about 5%. (I am embarrassed to admit I literally timing these things with my wristwatch.)

    Not sure if I understand your first point. I think I actually did what you suggested in the second version above.

  7. Pug says:

    I’ve started looking at Scala as a replacement for Matlab as well, so I have to say: this is pretty awesome. Thanks for posting it!

  8. rob says:

    Justin,

    I too have been trying to rid my research pipeline of matlab. The problem I find is that it’s great for some quick prototyping, but as soon as I want to do something that’s not a matrix operation, it quickly becomes a pain. I’d like to suggest another solution which I’ve been exploring, and with which I’ve found success… Python with Numpy/Scipy. Python is a very general and incredibly powerful language. It’s garbage collected, and allows for object oriented abstraction as well as functional programming tools. Perhaps one of the most interesting things about Python, and what makes it so tempting for numerical computing, is the ease of wrapping C/C++ and Fortran code.

    Numpy is a very comprehensive matrix library for Python. It allows most, if not all, of the operations afforded by Matlab, and uses a very natural syntax. Behind the scenes, it wraps C code for speed and efficiency. Scipy is a medley of other useful libraries and tools that deal with Numpy arrays as the basic array/matrix type. Together (and with matplotlib for plotting), they recreate almost all of the functionality of Matlab and many of the Matlab toolboxes. Furthermore, the numerical routines are fast as many scipy functions simply wrap calls to highly optimized fortran libraries. Finally, everything is exposed in Python with all the syntactic sugar and functional generality of the Python language. Anyway, if you’re still open to new options as Matlab replacements, I strongly suggest you take a look at Numpy and Scipy… I’ve found they do the trick wonderfully.

  9. justindomke says:

    Rob-

    I looked at python+numpy+scipy in great depth a few months ago. My impression was that it made an excellent matlab replacement. It had almost all the advantages of matlab plus actually being a nice language. The reason I ultimately decided not to switch is that if you can’t write the code in a vectorized form, python is slow. As in painfully, unacceptably slow. I’ll happily sacrifice a factor or 2 or whatever for more programmer efficiency, but the great language benchmarks game ( http://shootout.alioth.debian.org/ ) and my experience both suggest that the penalty is 100x-1000x. The reason I’m trying to ditch matlab is largely because I need fast, non-vectorizable, operations.

    I don’t mean this to criticize Python– the goals of the Python were never hardcore number crunching, and it doesn’t seem like bare-metal efficiency has been much of an issue to the user base. It would be very cool if the numpy/scipy community helped provided more incentive for projects like psyco or pypy or weave, because aside from speed, python wins almost every count. (But speed is a dealbreaker.)

  10. Zemian Deng says:

    Justin,

    Scala allow you do args(0).toInt instead of writing Integer.parseInt(args(0)).

    BTW, what tool you used for converting from pgm to png?

    I tried ToyViewer on Mac and it works good. http://www7a.biglobe.ne.jp/~ogihara/software/OSX/toyv-eng.html

    Thanks for the post!

  11. justindomke says:

    Zemian-

    I used IrfanView http://www.irfanview.com/ which, unfortunately, appears to be windows only.

  12. jvbsoares says:

    I just found Justin’s blog, and I’m really happy to see these suggestions!
    I’ve used Matlab a lot, for advantages that everybody already knows.
    My problems with it are:
    – not object-oriented (apparently they’re going to implement this, so this might not be a good excuse for not using it?)
    – it would be nice to have fast for-loops🙂
    – You guys in the US might find this strange, but I studied in a large Brazilian university and it was really a pain to get the license every month at our department; I’m not blaming this on Matlab, but there were times that we just couldn’t run it at our dept, and it was always a hassle; this is my major objection. And then everybody gets locked in, because everybody else has their code in Matlab; not good, IMHO.
    – I implemented a GUI in Matlab 6 and it hurt; a lot… It just wasn’t ready for developing GUIs. I heard Matlab 7 is a lot better for this, though…

    I’m going to try to play with Scala… and maybe do something directly in Java, just because I still like the fact that it’s ready to run on different platforms with no hassles; might not be worth it, though, I don’t know.

    best regards,
    Joao.

  13. jvbsoares says:

    Ok, so I just started reading a bit about Scala. It appears you can actually run a Scala program using the Java interpreter:
    http://www.codecommit.com/blog/scala/scala-for-java-refugees-part-1
    I’m pretty much sold…

    If I understood everything, Justin, it seems that you’re interested in the syntactic capabilites of Scala over Java? Are there other things I’m missing? I just started reading about this, so any pointers are welcome.

    thanks!
    Joao.

  14. jvbsoares says:

    http://www.javaworld.com/podcasts/jtech/2007/120607jtech007.html
    Bill Venners gives an instructional interview; answered all my questions, I think.

    regards,
    Joao.

  15. justindomke says:

    Hey Joao Honestly, I would probably just work in Java, except for the fact that it has no operator overloading. I refuse to write “A.set(2,5,mult(B.get(1,2),C.get(2,7)))” instead of “A(2,5)=B(1,2)*C(2,7)”, which is what Java would seem to demand.

  16. John Kimball says:

    Justin;

    In your explorations of python/numpy/scipy did you look at sage and cython?

    Sage (http://www.sagemath.org/) integrates python, numpy, scipy and many other open source packages into an integrated python-based environment intended to provide a viable alternative to Matlab and the like.

    Regarding the slow non-vectorizable Python code, the Sage project is developing a fork of the pyrex project called Cython (http://www.cython.org) that allows for loops to be executed in C code with very little programmer effort. Cython can be combined with numpy. I’ve seen claims of speeds approaching C/C++, though I haven’t used Cython yet.

    That said, I’ll be following your explorations of numerical computation in Scala with interest.

    John

  17. justindomke says:

    John-

    Yes, I seriously considered the sage environment. If I wasn’t worried about efficiency, that’s probably what I would use. On paper, cython seems like a fairly nice solution to the speed problem. However, when I tried it out, I was never able to get cython to run at anything near C speed for anything other than trivial examples. It is probable that with the appropriate tricks it would be as fast, but there wasn’t nearly enough documentation for me to figure out what I was doing wrong, so I gave up. (It would be easier to just code in pure C, where everything runs at exactly the speed you expect it to.)

    Another option is python+weave, which allows you to inline c code *in the same function*, with almost no boilerplate code. Still, this runs into the same basic problem as my current setup (matlab + c for the slow stuff), namely that sometimes it just isn’t possible to make a clean code separation between high level abstract code, and low level code that needs to be run at high speed. If you can’t cleanly interface between the two, any solution based on two languages (pyton+cython, python+weave, matlab+c, etc.) falls apart.

    Thanks for your comment,
    Justin

  18. James says:

    Hi Justin,

    Do you have any speed comparisons for your Mandelbrot code? Does this run much slower than Java or C?

  19. Tim says:

    Have you seen this livejournal.
    http://ignaciopaulyx.livejournal.com/2009/04/05/
    Looks like you have a fan/plaigarist.

  20. dramage says:

    I’m the author of Scalala and it’s kind of gotten to the point where it’s reasonably stable, so I’ve written up a quick start guide linked from the scalala googlecode page.

  21. Ewald says:

    Hi,

    Not sure if this is relevant anymore. Just for fun I coded your second example in python and cython (the calculation bit).

    Using your numbers, n=2048 level=100 I get:

    Python Mandelbrot with n=2048 and level=100 in 58.295 seconds
    Cython Mandelbrot with n=2048 and level=100 in 1.218 seconds
    Same answer? : True
    Speedup: 47.8744235634

    Laptop has Core 2 Duo T9300 with 3GB ram on Ubuntu 8.10 64bit. I can send you the code if you’d like. Using python/cython/scipy/numpy/matplotlib/f2py is probably a better way than more obscure (for number crunching and real world problems) languages.

    Cheers

  22. Justin says:

    Very cool. Actually, since writing this post, I have basically abandoned Scala for python, so that is encouraging. (Python’s libraries are 1000x better, which outweighs everything in the end) I am frankly amazed you got such good performance out of regular python! I’ve been using weave, rather than cython for performance, but I should probably take another look at that. I would appreciate looking at the code ((myfirstname).(mylastname).rit.edu, or just post it as a comment). I remember having a hard time figuring out how to get cython to run quickly.

  23. Ewald says:

    http://wiki.cython.org/tutorials/numpy

    Look here for info on how to use numpy arrays effciently inside a cython routine. I mailed the code.

    cython -a yourfile.pyx makes an html file that highlights the slow (python calls) bits. Double click the lines to see the generated code.

  24. Brett England says:

    The while loop can be removed if you use a tail recursive function. This is the functional approach that Scala aspires you to use:

    def iterate(z:Complex, c:Complex, level:Int, i:Int): (Complex,Int) =
    if(z.abs > 2 || i > level) (z,i) else iterate(z*z+c, c, level, i+1)


    val (z, i) = iterate(new Complex(0,0), new Complex(x,y), level, 0)
    out.write(if (z.abs < 2) 0 else i)

  25. bmdesignhki says:

    Interesting. I’m in 2013 and the same situation seems to go on. Scala would be a splendid language in many ways (as described above) but there’s no obvious plotting library (anyone please correct me if I’m wrong! Scalala in particular looks dead-ish.).

    Eventually, someone will make one out of JavaFX (for example). It can be used not only to make pie chart GUIs. In the mean time, anything?

    I’ve moved 100% to Scala coding last summer. Now I end up learning matplotlib to make a couple fancy pictures. If *anyone* starts making a plotting library for Scala, count me in as a user / tester / motivator.

  26. uday says:

    great post. I’ve been looking for a Matlab replacement and after trying out a lot of languages settled down on Scala. i have a serious requirement for easy parallel programming and Matlab’s efficiency of simply changing “for” to “parfor” is almost unbeatable – only Scala has a similar functionality.

    R – too slow, and unbelievably too buggy, but lovely syntax (perhaps the best of all languages). tons of parallel programming packages, none of which work and the few that do end up spending more time than a non-parallel code.

    Python – great if you have a laundry list of small talks. can’t ever understand how others can one implement complicated algorithms. this weird one instruction on one line and the tabs creates too long files that I get lost in my own code. error detection is too painful, debugging takes foreover, very poor memory handling for large data (relative to even R!), and of course no backward compatibility, python 2.5 codes don’t run in 2.6 which don’t run in 2.7 and now python 3.3 is an entirely different syntax which also btw doesn’t run 3.2 or 3.1 codes.

  27. Python’s universal indentation makes for great readability, and PyPy/ShedSkin/Cython/CorePy can exceed C speed, but i hate its mandatory self parameters in classes. Scala, using the fast JVM, also solved that language wart.

    @uday, Python 2.7 runs most Python 2 and some Python 3 code. Python 3 only runs Python 3 code; this mostly means you’ll need to add () to print and change Unicode handling, which should be automated by the 2to3 tool: http://docs.python.org/3/library/2to3.html

  28. Pingback: Mandelbrot集合(多图超能预警)格物致知 | 格物致知

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s