Isolated Mandelbrot Set Explorer

mandelbrot-equation

It amazes me that an equation so simple can produce images so beautiful and complex.  I was first introduced to the Mandelbrot set by the August 1985 issue of Scientific American. Wikipedia says “The Mandelbrot set is the set of complex numbers ‘c’ for which the sequence ( c , c² + c , (c²+c)² + c , ((c²+c)²+c)² + c , (((c²+c)²+c)²+c)² + c , …) does not approach infinity.” This means that in the images below, the areas in black are points in the Mandelbrot set and the other point colors are determined by how quickly that point “zooms” off towards infinity. The links above give a deeper discussion on the Mandelbrot set – for this post I’d like to focus on how I was able to speed my Mandelbrot set explorer up using isolates.

Mandelbrot Set
I find I learn things more easily through application than abstraction. In 2006 I wanted to learn more about the built-in GUI features of Dyalog for Windows and the application I chose was a Mandelbrot set explorer. While the basic algorithm for calculating a point in the Mandelbrot set is fairly simple, a 640×480 pixel image could involve over 75 million calculations. As new versions of Dyalog are released, I try to experiment with any applicable new features and possibly incorporate them into the explorer. For instance, when complex numbers were introduced in Dyalog, I modified the explorer to use them.

Dyalog v14.0 introduced a prototype of a parallel computing feature called isolates that may one day be integrated into the interpreter. Isolates behave much like namespaces except that they run in their own process with their own copy of the interpreter. Each process can run on a separate processor, greatly improving performance. You can even set up isolate servers on other machines.

Mandelbrot set calculations are ideal candidates for parallelization as:

  • they’re CPU intensive
  • each pixel can be calculated independently of others
  • there are no side effects and no state to maintain

The core calculation is shown below. It takes the real and imaginary parameters for the top left corner, the height and width of the image and the increment along each dimension and uses those parameters to build a set of co-ordinates, then iterates to see how quickly each co-ordinate “escapes” the Mandelbrot set.

 ∇ r←real buildset_core imaginary;height;i_incr;top;width;r_incr;left;set;coord;inds;cnt;escaped
  ⍝ Calculate new Mandelbrot set - remember using origin 0 (⎕IO←0)
  ⍝ real and imaginary are the parameters of the space to calculate
  ⍝       real            imaginary
  ⍝ [0]   width           height            of the image
  ⍝ [1]   increment size  increment size    along each axis
  ⍝ [2]   top             left              coordinates

   (height i_incr top)←imaginary            ⍝ split imaginary argument parameters
   (width r_incr left)←real                 ⍝ split real argument parameters
   set←coord←,(0J1×top+i_incr×⍳height)∘.+(left+r_incr×⍳width) ⍝ create all points in the image
   inds←⍳≢r←(≢coord)⍴0                      ⍝ initialize the result and create vector of indices
   :For cnt :In 1↓⍳256                      ⍝ loop up to 255 times (the size of our color palette)
       escaped←4<coord×+coord               ⍝ mark those that have escaped the Mandelbrot set
       r[escaped/inds]←cnt                  ⍝ set their index in the color palette
       (inds coord)←(~escaped)∘/¨inds coord ⍝ keep those that have not yet escaped
       :If 0∊⍴inds ⋄ :Leave ⋄ :EndIf        ⍝ anything left to do?
       coord←set[inds]+×⍨coord              ⍝ the core Mandelbrot computation... z←(z*2)+c
   :EndFor
   r←(height,width)⍴r                       ⍝ reshape to image size
 ∇

So, what’s involved in using isolates to parallelize the explorer?  As it turns out, not much! First I )COPYed the isolate namespace from the isolate workspace distributed with Dyalog v14.0. I wanted to be able to select the number of processors to use, so I added a line to my init function to query how many processors are available:

 processors←1111⌶⍬ ⍝ query number of processors

Next I wrote a small function to initialize isolates for as many processors as were selected (2 or more). Notice that when creating the isolates, I copy buildset_core (from above) into each one:

 ∇ init_isolates processors
  ⍝ called upon initialization, or whenever user changes number of processors to use
   f.isomsg.Visible←1    ⍝ turn form's "Isolates initializing" message on
   {}isolate.Reset''     ⍝ reset the isolates
   {}isolate.Config'processors'processors ⍝ set the number of processors
   :If processors>1      ⍝ if using more than one processor
       ⍝ create the isolates, populating each one with buildset_core
       isolates←isolate.New¨processors⍴⎕NS'buildset_core'
   :EndIf
   f.isomsg.Visible←0    ⍝ turn the form message off
 ∇

Finally, I wrote a routine to share the work between the selected number of processors and assemble the results back together.

 ∇ r←size buildset rect;top;bottom;left;right;width;height;vert;horz;rows;incr;imaginary;real
  ⍝ Build new Mandelbrot Set using isolates if specified
  ⍝ rect is the top, bottom, left, and right
  ⍝ size is the height and width of the image
   (top bottom left right)←rect
   (height width)←size
   vert←bottom-top                             ⍝ vertical dimension
   horz←right-left                             ⍝ horizontal dimension
   rows←processors{¯2-/⌈⍵,⍨(⍳⍺)×⍵÷⍺}height     ⍝ divide image among processors (height=+/rows)
   incr←vert÷height                            ⍝ vertical increment
   imaginary←↓rows,incr,⍪top+incrׯ1↓0,+\rows  ⍝ imaginary arguments to buildset_core
   real←width,(horz÷width),left                ⍝ real arguments to buildset_core
   :If processors=1                            ⍝ if not using isolates
       r←real buildset_core⊃imaginary          ⍝ build the image here
   :Else
       r←⊃⍪/(⊂real)isolates.buildset_core imaginary ⍝ otherwise let the isolates to it
   :EndIf
 ∇

The images below depicts what happens with 4 processors. Each processor builds a slice of the image and the main task assembles them together.
Drawing1

So, how much faster is it?  Did my 8 processor Intel Core i7 speed things up by a factor of 8? No, but depending on where in the Mandelbrot set you explore, the speed is generally increased by at least a factor of 2 and sometimes as high as 5 or 6. I thought about why this might be and have a couple of ideas…

  • There’s overhead to receive and assemble the results from each isolate.
  • The overall time will be the maximum of the individual times; if one slice is in an area that requires more iterations, then the main task needs to wait for those results, even if the other slices have already finished. I could possibly show this by updating the individual slices as they become available. Maybe next time 🙂

What impressed me was that I could speed up my application significantly with about 20 lines of code and maybe an hour’s effort – most of which was spent tinkering to learn isolates. You can download the Mandelbrot set explorer from my GitHub repository.

Making Controlled Turns with the DyaBot

This blog originally started when I took delivery of the DyaBot, a Raspberry Pi and Arduino based variant of the C3Pi running Dyalog v13.2. The architecture of the ‘bot and instructions for building your own inexpensive robot can all be found in blog entries from April to July of last year.

The downside of only using inexpensive components is that some of them are not very precise. The worst problem we face is that the amount of wheel movement generated by the application of a particular power level varies from one motor to the next, and indeed from moment to moment. Trying to drive a specific distance in a straight line, or make an exact 90 degree turn regardless of the surface that the ‘bot is standing on, are impossible tasks with the original Dyabot. You can have a lot of fun, as I hope the early posts demonstrate, but we have higher ambitions!

Our next task is to add motion sensors to the DyaBot, with a goal of being able to measure actual motion with sufficient accuracy to maintain our heading while driving straight ahead – and to make exact turns to new headings, like the 90-degree turn made in this video (the ‘bot has been placed on a Persian carpet to provide a background containing right angles):

Introducing the MPU-6050

For some time we have had an MPU-6050, which has 3-axis rotation and acceleration sensors, attached to our I2C bus, but haven’t been able to use it. About ten days ago (sorry, I took a few days off last week), @RomillyC came to visit us in Bramley to help me read some Python code that was written for this device. The current acceleration or rotation on each axis is constantly available as a register and can be queried via I2C. Translated into APL, the code to read the current rate of rotation around the vertical axis is:

 ∇ r←z_gyro_reading ⍝ current z-axis rotation in degrees
   r←(read_word_2c 70)÷131 ⍝ Read register #70 and scale to deg/sec
 ∇

 ∇ r←read_word_2c n;zz
   r←I2 256⊥read_byte_2c¨n+0 1               ⍝ Read 2 registers and make signed I2
 ∇

 ∇ r←read_byte_2c n;z
   z←#.I2C.WriteArray i2c_address (1,n)0     ⍝ Write register address
   r←2 2⊃#.I2C.ReadArray i2c_address (1 0)1  ⍝ Read input
 ∇

 I2←{⍵>32767:1-65535-⍵ ⋄ ⍵} ⍝ 16-bit 2's complement to signed

Integrating to Track Attitude

Reading the register directly gives us the instantaneous rate of rotation. If we want to track the attitude, we need to integrate the rates over time. The next layer of functions allows us to reset the attitude and perform very primitive integration by multiplying the current rotation with the elapsed time since the last measurement:

∇ r←z_gyro_reset            ⍝ Reset Gyro Integration
  z_gyro_time←now           ⍝ current time in ms
  z_gyro_posn←0             ⍝ current rotation is 0
∇

∇ r←z_gyro;t;delta          ⍝ Integrated gyro rotation
  t←now                     ⍝ current time in ms
  delta←0.001×t-z_gyro_time ⍝ elapsed time in seconds
  r←z_gyro_posn←z_gyro_posn+delta×z_gyro_reading+z_gyro_offset
  z_gyro_time←t             ⍝ last time
∇

We can now write a function to rotate through a given number of degrees (assuming that the gyro has just been reset): we set the wheels in motion, monitor the integrated angle until we are nearly there, then shut off the engines. We continue logging the angles for a brief moment to monitor the deceleration. The function returns a two-column matrix containing the “log” of rotation angles and timestamps:

    ∇ log←rotate angle;bot;i;t
     
      log←0 2⍴0 ⍝ Angle, time
      bot.Speed←3 ¯3 ⍝ Slow counter-clockwise rotation
     
      :Trap 1000 ⍝ Catch interrupts
          :Repeat
              log⍪←(t←z_gyro),now
          :Until t>angle-7 ⍝ Stop rotating 7 degrees before we are done    
      :EndTrap

      bot.Speed←0 ⍝ cut power to the motors
      :For i :In ⍳25 ⍝ Capture data as we slow down
          log⍪←z_gyro,now
      :EndFor
    ∇

The logged data from the rotation captured in the video at the top is charted below:

Controlled 90-Degree RotationEven with the primitive integration and rotation strategies, the results already look quite promising. I’ll be taking most of this week off – part of it without access to the internet(!), but once I am back, expect the next blog entry to explore writing functions that accelerate and slow down in a more controlled fashion as well as stop at the right spot rather than relying on a specific amount of friction to rotate through the last 7 degrees (note the very slight reverse rotation at the end, probably caused by the Persian carpet being a bit “springy”). I will also clean up the code and post the complete solution on GitHub – and perhaps even look at some better integration techniques.

If you would like to make sure you don’t miss the next installment, follow the RSS feed or Dyalog on Twitter.

Dancing with the Bots

Last week the ‘bots were busy preparing for the J Language Conference in Toronto, where they made their first public appearance together. Upon returning to Bramley they continued training and we are proud to present the first recording of their new dance:

The ‘bots are both running the same DyaBot class as last year. This class exposes a property called Speed, which is a 2-element vector representing the speed of the right and left wheels respectively. Valid values range from +100 (full speed ahead) to -100 (full reverse). The annotations displayed at the top left show the settings used for each step of the dance.

Controlling Two Robots at Once using Isolates

Isolates are a new feature included with Dyalog version 14.0, designed to make it easy to perform distributed processing. In addition to making it easy to used all the cores on your own laptop or workstation, isolates make it possible to harness the power of other machines. This requires the launching of an “isolate server” on each machine that wants to offer its services:

Starting an isolate server using PuTTY to run Dyalog on the robot.

Starting an isolate server on DyaBot00 using PuTTY.

Once we have an isolate server running on each robot we can take control of them from a remote session as follows:

      )load isolate
      #.isolate.AddServer 'dyabot00' (7052)      
      #.isolate.AddServer 'dyabot04' (7052)
      bots←isolate.New¨Bot Bot
      bots.Init
 dyabot00  dyabot04

Above, we create two instances of the Bot namespace. The expression Bots.Init invokes the Init function, which returns the hostname, in each isolate:

:Namespace Bot

    ∇ r←Init;pwd
      pwd←∊⎕SH'pwd' ⍝ Find out where to copy from
      #.⎕CY botws←pwd,'/DyaBot/DyaBot.dws' ⍝ copy ws
      i←⎕NEW #.DyaBot ⍬ ⍝ Make DyaBot instance    
      r←⎕SH'hostname' ⍝ Return hostname
    ∇

:EndNamespace

Next, we define a function “run” that will take a vector of dance steps as input. Each step is a character vector (because that makes editing slightly easier!) containing five numbers: The first two set the speed of one robot, the next two the speed of the other and the fifth defines the duration of the step. After each step we pause for a second, to give humans time to appreciate the spectacle:


    ∇ run cmds;data;i;cmd;z
[1]    ⎕DL 5
[2]    :For i :In ⍳≢cmds
[3]        :If ' '∨.≠cmd←i⊃cmds
[4]            data←1 0 1 0 1⊂2⊃⎕VFI cmd ⍝ Cut into 3 numeric pieces
[5]            z←bots.{i.Speed←⍵}2↑data ⋄ ⎕DL⊃¯1↑data ⋄ z←bots.(i.Speed←0)
[6]            ⎕DL 1
[7]        :EndIf
[8]    :EndFor
    ∇

Now we are ready to roll: Call the run function with a suitable array and watch the robots dance (see the video at the top)!

      ↑choreography
50  50  0   0 1.5
 0   0 50  60 1.2
50 ¯50 50 ¯50 0.3
20  80 10  70 5  
50 ¯50 50 ¯50 0.3
50  50  0   0 1.5
 0   0 50  60 1.2

      dance choreography

Join us again next week to hear what happened when Romilly came to Bramley to help wire up the accelerometer and gyro!