Parallel mandelbrot generator

F# offers very good support for implementing algorithms in a parallel way to exploit multi-core processors. The following example of a mandelbrot generator uses a lot of F# specific features

  1. using operator overloading to create a complex type,
  2. the async keyword to implement parallel computation,
  3. the forward pipe operator ‘|>’ to join multiple steps of computations, and
  4. the System.Drawing package to create the output image.

Lets first create a complex type to allow operation in the complex number space.

> type complex(r:float, i:float) =
-   member this.r = r
-   member this.i = i
-
-   static member(+) (x:complex, y:complex) =
-     complex(x.r+y.r, x.i+y.i)
-
-   static member(*) (x:complex, y:complex) =
-     complex(x.r*y.r-x.i*y.i, x.r*y.i+x.i*y.r)
-
-   static member r2 (x:complex) = x.r*x.r + x.i*x.i
-
- ;;

type complex =
 class
 new : r:float * i:float -> complex
 member i : float
 member r : float
 static member ( + ) : x:complex * y:complex -> complex
 static member ( * ) : x:complex * y:complex -> complex
 static member r2 : x:complex -> float
 end
> complex(1.0, 2.0) * complex(3.0, 4.0);;
val it : complex = FSI_0240+complex {i = 10.0;
 r = -5.0;}
>

Next comes the mandelbrot algorithm. This algorithm takes a complex number and the maximum number of iterations to compute if this complex number is part of the mandelbrot set (=1) or not (=0). The implementation is a tail-call recursive implementation using the embedded function ‘iterate’

> let mandelbrot c niter =
-   let rec iterate z i =
-     let zn = z * z + c
-     if complex.r2(zn) > 2.0 then 0
-     else
-       if i = 0 then 1
-       else iterate zn (i-1)
-   iterate (complex(0., 0.)) niter
-
- ;;

val mandelbrot : complex -> int -> int

> mandelbrot (complex(0.3, 0.5)) 50;;
val it : int = 1

The function seq2bitmaprow is a helper function to write a sequence (e.g. [|0, 1, 1, 1, 0|]) into a bitmap row. The sequence is assumed to have the same number of elements as the bitmap width in pixels.

It is worth mentioning that the last argument of seq2bitmaprow is the sequence, this will allow us later to pipe the result sequence into this function.

> let seq2bitmaprow (bitmap:Bitmap) row seq =
-   let mutable i = 0
-   for x in seq do
-     if x = 1 then
-       bitmap.SetPixel(i, row, Color.Black)
-     else
-       bitmap.SetPixel(i, row, Color.White)
-     i <- i + 1
-
- ;;

val seq2bitmaprow : Bitmap -> int -> seq<int> -> unit

The next part is about plumbing it all together using the forward pipe operator. Try out the following code snippet in the interactive to see how the row sequence [|0..(w-1)|] is parallely piped through the mandelbrot algorithm.

> let w = 50
- let h = 50
- let y = 25;;

val w : int = 50
val h : int = 50
val y : int = 25

> let row = [|0..(w-1)|];;

val row : int [] =
 [|0; 1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16; 17; 18; 19; 20;
 21; 22; 23; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39;
 40; 41; 42; 43; 44; 45; 46; 47; 48; 49|]

> let async_comp = row |> Array.map
-                         (fun x -> async {
-                           return (mandelbrot (complex(float (x*2)/(float w)-1.5,
-                                                       float (y*2)/(float h)-1.0)) 50- )})
- ;;

val async_comp : Async<int> [] =
 [|Microsoft.FSharp.Control.FSharpAsync`1[System.Int32];
 Microsoft.FSharp.Control.FSharpAsync`1[System.Int32];
 Microsoft.FSharp.Control.FSharpAsync`1[System.Int32];
...

> let para_comp = async_comp |> Async.Parallel;;

val para_comp : Async<int []>

> let row_result = para_comp |> Async.RunSynchronously;;

val row_result : int [] =
 [|0; 0; 0; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1;
 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 0; 0; 0; 0; 0; 0|]

Full source:

open System
open System.IO
open System.Drawing

type complex(r:float, i:float) =
  member this.r = r
  member this.i = i

  static member(+) (x:complex, y:complex) =
    complex(x.r+y.r, x.i+y.i)

  static member(*) (x:complex, y:complex) =
    complex(x.r*y.r-x.i*y.i, x.r*y.i+x.i*y.r)

  static member r2 (x:complex) = x.r*x.r + x.i*x.i

let mandelbrot c niter =
  let rec iterate z i =
    let zn = z * z + c
    if complex.r2(zn) > 2.0 then 0
    else
      if i = 0 then 1
      else iterate zn (i-1)
  iterate (complex(0., 0.)) niter

let seq2bitmaprow (bitmap:Bitmap) row seq =
 let mutable i = 0
 for x in seq do
   if x = 1 then
     bitmap.SetPixel(i, row, Color.Black)
   else
     bitmap.SetPixel(i, row, Color.White)
   i <- i + 1

let () =
 let w = int(Environment.GetCommandLineArgs().[1])
 let h = w
 let bitmap = new Bitmap(w, h)
 let row = [|0..(w-1)|]
 for y in [0..(h-1)] do
   [|0..(w-1)|]
     |> Array.map
       (fun x -> async {
         return (mandelbrot (complex(float (x*2)/(float w)-1.5,
                                     float (y*2)/(float h)-1.0)) 50)})
     |> Async.Parallel
     |> Async.RunSynchronously
     |> seq2bitmaprow bitmap y

 bitmap.Save("mandel.png")

Resulting image (100×100):

Further reading: Asynchronous Workflows (F#)

Advertisements
This entry was posted in Algorithms, Parallel. Bookmark the permalink.

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