Hacking with APL

Vassar Hackathon Poster

Vassar Hackathon Poster

Thanks to our dear friend Dr. Ray Polivka, Dan Baronet and I had the opportunity to travel to Vassar College to participate in their Community Hackathon held on 5-6 February 2016.

“What’s a hackathon?” you ask?
Well, we did too, as we’d never participated in one before.  :)
According to the Hackathon’s announcement:

“CommunityHack is a way to bridge the gap between Vassar CS Majors and Vassar students in other departments as well as local high school and community college students in order to provide them with the opportunity to explore innovative applications of Computer Science. Dance, music, art, video games? CS can be incorporated into all of that and more!”

StuCommunityHack_Sponsorsdents from Vassar as well as nearby colleges and high schools were invited to attend. In other words, it was a great opportunity to introduce APL to a new generation.  As this was our first Hackathon, we had no idea what to expect.  Laura, the Hackathon’s organizer, did a wonderful job organizing the event and making us feel welcome. We were invited to give an hour long presentation and Dyalog was listed as an event sponsor.

The Hackathon was a 24 hour event where students were encouraged to split up into groups and pick a problem to solve.  During the course of the event, presentations were made on a variety of subjects including “Autonomous Robots Test Ideas About the Evolution of Brains”, “How to make games quick!”, “Virtual Reality”, and of course “Hacking with APL”. Friday evening started with introductions and ice-breakers. During our introduction, I was able to talk a bit about APL and the presentation we would be making on Saturday. Apparently this generated some interest as a number of students came up to Dan, Ray, Jon McGrew and me to ask about APL. We spent several hours showing them APL, to which they seemed eagerly receptive.

I had the pleasure of working with Grace, a CS sophomore, to implement APL (A Puppy Language) in APL. Her project idea was to write an application for potential puppy owners to use so they could get an idea of the responsibility of owning and caring for a puppy. We worked into the wee hours of the night and wound up implementing a multi-threaded domain-specific-language (DSL) where the “puppy”, running in a separate thread, would react to commands typed into the APL session. Negative actions and ignoring the puppy would cause the the puppy’s happiness points (PHPs) to decrease whereas positive actions would increase PHPs.   Grace seemed to really enjoy working with APL and returned to the hackathon twice on Saturday as her schedule permitted to continue work on her project.

Saturday, I was slightly concerned that following a talk on virtual reality, APL might not seem all that “cool”, but my fears were allayed for, as I was waiting before my presentation, several students asked the person at the registration desk specifically about the APL presentation.

HackingWithAPL

The presentation went rather well.  Watching the jaw-dropping “Wow!” expressions on the faces of many of the students as we showed the power of APL notation made me reminisce back to the time when I first learned APL and was amazed at what I could do with a computer.  It also reminded me how blessed I’ve been to have used APL throughout my career.

Our participation in the Hackathon was a great experience. We were able to show Dyalog to close to 100 students, promote the upcoming APL Problem Solving Competition, and encourage people to download and try Dyalog – we had 18 student downloads over the Hackathon weekend. This may have been our first Hackathon, but I’m certain it won’t be our last.

On a personal note, after Dan and I drove up to Montreal to spend the upcoming week working with the APL Tools Team, I received a very nice email from Grace where she wrote “I just wanted to thank you so much for taking the time to work with me on puppy.dws – it is currently my favorite thing that I have ever made.” and “It was really fun working in APL, and I will definitely check out the Dyalog competition.”

HackingWithAPL3

 

Solving the 2014 APL Problem Solving Competition – How Tweet It Is

This post is the continuation of the series where we examine some of the problems selected for the 2014 APL Problem Solving Competition.

The problems presented in Phase 1 of the competition were selected because they could be solved succinctly, generally in a single line of APL code. This makes them well suited for experimentation on TryAPL.org.

Problem 2 of Phase 1, entitled “How tweet it is” reads

“Twitter messages have a 140 character limit – what if the limit was even shorter? One way to shorten the message yet retain most readability is to remove interior vowels from its words. Write a dfn which takes a character vector and removes the interior vowels from each word.”

Test cases:
      {your_solution} 'if you can read this, it worked!'
if yu cn rd ths, it wrkd!
      {your_solution} 'APL is REALLY cool'
APL is RLLY cl
      {your_solution} '' ⍝ an empty vector argument should return an empty vector

      {your_solution} 'a' ⍝ your solution should work with a single character message
a

We’ll examine a couple of approaches to this problem – one that’s more “traditional APL” code, and another that makes use of a really helpful Dyalog feature.

This problem could be restated as “find and remove the vowels that aren’t at the beginning or end of a word”. To start with, we need to determine where the words are and where the vowels are. A word is a contiguous set of letters; multiple words are separated by spaces or punctuation. For simplicity’s sake, we’ll ignore contractions and possessives.

The “Traditional APL” Approach

This approach employs a technique that is not commonly found outside of APL and its brethren – using a Boolean vector to determine which elements to remove or keep. First, let’s find where all the vowels are:

      string←'If you can read this, it worked!'
      vowels←{⍵∊'aeiouAEIOU'}
      vowels string
1 0 0 0 1 1 0 0 1 0 0 0 1 1 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 0

To help illustrate what’s happening, I’ll write a little operator called “show” to compactly display the string, the Boolean vector, and the elements that would be selected by applying the Boolean to the string.

      show←{⍵⍪⍵{↑(1 0⍕ ⍵)(⍵\⍵/⍺)}⍺⍺ ⍵}
      vowels show string
If you can read this, it worked!
10001100100011000010001000100100
I   ou  a   ea    i   i   o  e

Next we want to remove vowels that aren’t at either end of a word. First, find where the words are by finding where the letters are.  There are several ways to do this; the most obvious may be to use a character vector constant:

      letters←{⍵∊'abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ'}

Long character constants seem a bit awkward to me.  So, another technique uses the Unicode Conversion system function to return the 26 characters starting at the code points for each of ‘a’ and ‘A’:

      letters←{⍵∊⎕UCS (⍳26)∘.+¯1+⎕UCS'aA'}

Yet another way might be to use the code point values directly and do numerical operations:

      letters←{{((⍵≥65)∧⍵≤90)∨(⍵≥97)∧⍵≤122}⎕UCS ⍵}

Which technique you choose is largely a matter of taste and style. All three return the same result and have comparable performance. My personal preference is the second one – it has fewer characters for me to mistype :)

      letters show string
If you can read this, it worked!
11011101110111101111001101111110
If you can read this  it worked 

So now let’s mark the interior letters of the words. This employs a technique known as shift and compare that I learned in the early 1980s when I was privileged to work with Bob Smith. Among Bob’s many contributions to the APL world was a book on Boolean Functions and Techniques. To mark the interior letters, we’ll do both a right and left shift:

      interior←{⍵∧(¯1↓0,⍵)∧1↓⍵,0}
      {interior letters ⍵} show string
If you can read this, it worked!
00001000100011000110000000111100
    o   a   ea   hi       orke  

The last step is to find interior vowels and negate:

      {(vowels ⍵)∧interior letters ⍵} show string
If you can read this, it worked!
00001000100011000010000000100100
    o   a   ea    i       o  e  

      {(vowels ⍵)⍲interior letters ⍵} show string
If you can read this, it worked!
11110111011100111101111111011011
If y u c n r  d th s, it w rk d!

Putting it all together…

      tweet←{⍵/⍨~(⍵∊'aeiouAEIOU')∧{(1↓⍵,0)∧¯1↓0,⍵}⍵∊⎕UCS(⍳26)∘.+¯1+⎕UCS'aA'}
      tweet string
If yu cn rd ths, it wrkd!

The Other Technique – Using Regular Expressions

In version 13.0, Dyalog introduced the system functions ⎕S and ⎕R as interfaces to the PCRE (Perl Compatible Regular Expression) library. Like APL, regular expressions may seem a bit alien at first, but in the years since their incorporation into Dyalog, I’ve grown to appreciate their power and flexibility – they can frequently accomplish complex string manipulations more succinctly than their APL equivalents thus furthering Dyalog’s power as a tool of thought, notation and execution.

      tweet←{('\B[AEIOU]\B' ⎕R '' ⍠ 1) ⍵}
      tweet string
If yu cn rd ths, it wrkd!

The expression above replaces any vowel (⍠ 1means case-insensitive) that is not at the beginning or end of a word with the empty vector, effectively removing the interior vowels. A blog post is not enough space to give an adequate overview of regular expressions. But I hope the expression above piques your interest and encourages you to experiment with ⎕S and ⎕R on TryAPL.org or with a Dyalog system of your own.

Simply A-maze-ing

maze

One of many things I like about APL is that it’s fun to use for recreational computing. I will frequently happen upon an interesting problem, puzzle, or piece of code and consider how I might implement it in APL.

I was thinking about how to generate mazes for an idea I have about a game to help kids learn APL (that may be a topic for a future blog entry). Anyhow, I found an interesting web page about maze generating algorithms where the author found one, the “Growing Tree Algorithm“, to be of particular interest. His page included roughly 100 lines Ruby code to implement the algorithm. The algorithm can be boiled down to:

  • Seed a list of visited cells with a cell selected at random
  • While there are unvisited cells
    • If the current cell has any unvisited neighboring cells
      • Select one at random
      • Remove the wall between the cells
      • Add the new cell to the list of visited cells
    • Otherwise backtrack (drop from the visited cell list) until you find a cell with an unvisited neighbor

Here’s a clip of the algorithm implemented in APL building a 10×10 maze.
Notice how whenever it hits a “dead end” it backtracks to the last cell that hasn’t been visited.

What might an APL approach to this algorithm look like? How to represent the maze? My first thought was to separate the maze generation from the drawing. After an hour (or so ;-) ) of tinkering, I’d come up with something that seemed to work pretty well and took about a dozen lines of code.

When I originally thought about writing this blog entry, I was going to launch into a discussion of the code and I realized that it might get lengthy and (egad!) boring. So instead, I’ll highlight a couple of the clever bits, show you the core maze generation code, and point you at the complete namespace for your own amusement and experimentation.

First the clever bits (at least I hope they’re clever)…

  • Represent the cells of the maze as an integer matrix where each element is an encoding of the walls surrounding each cell.  Use powers of 2 for the encoding.
  • Precalculate the indices of the neighboring cells around each cell so the core loop only has to use indexing and no on-the-fly computation.
  • Write a function to remove the common wall between two cells. I originally named the function “Reagan” (after President Reagan’s 1987 exhortation “Mr. Gorbachev, tear down this wall”), but in the spirit of political mindfulness, I renamed it “dewall”.

The core code for the maze generation looks like this:

∇ cells←maze size;valid;neighbors;dewall;visited;current;n;choices;next
 ⍝ Maze - modeled from http://weblog.jamisbuck.org/2011/1/27/maze-generation-growing-tree-algorithm
 ⍝ BPB  - Dec2014
 ⍝ size  - size of the maze in rows/cols
 ⍝ cells   - (2⍴size) integer matrix describing the walls around each cell using powers of 2
 ⍝       1
 ⍝     ┌───┐         ┌───┐
 ⍝   8 │   │ 2       │   │ = 11 = 1 + 2 + 8
 ⍝     └───┘
 ⍝       4
  size←2⍴size  ⍝ set size to square if singleton supplied as argument     

  valid←{{⍵/⍨⍵≢¨⊂⍬}size∘{(0∊⍵)⍱⍵∨.>⍺:⍵ ⋄ ⍬}¨⍵} ⍝ determine if a maze coordinate is valid
  neighbors←valid¨↓(⍳size)∘.+,1 ¯1∘.×1 0⌽¨⊂1 0 ⍝ calculate neighbors for each cell
     
  dewall←{{⍵[2]=0:{(1=⍵)⌽4 1}⍵[1] ⋄ {(1=⍵)⌽2 8}⍵[2]}⍺-⍵}  ⍝ remove wall between cells
     
  cells←size⍴15 ⍝ all cells start with 4 walls
     
  visited←,⊂?size ⍝ random starting cell
     
  :While 15∊cells       ⍝ while we still have cells to examine
      current←1↑visited   ⍝ pop the most recent cell
      :If 0=n←⍴choices←cells{⍵/⍨15=⍺[⍵]}current⊃neighbors ⍝ does it have any unvisited neighbors?
          visited↓⍨←1       ⍝ if none, backtrack
      :Else
          visited,⍨←next←choices[?n] ⍝ otherwise, add the new cell to the front of the list
          cells[next current]-←⊃next⊃.dewall current ⍝ update cell values for which wall was removed
      :EndIf
  :EndWhile
∇

You can get all the code in my maze generating namespace from GitHub. Save a local copy and use the SALT Load command to load it into your workspace, or just cut and paste it into your own namespace script with the editor. The maze namespace contains the following functions of interest:

 

      intmat←{animate} maze size
  • animate is an optional Boolean to indicate whether to animate the maze generation
  • size is the size (rows,cols) of the maze to generate; a single number generates a square maze
  • intmat is the integer matrix representation of the maze

For example: mat←1 maze 10 animates and generates a 10×10 maze

 

      z←{line} draw intmat
  • line is an optional Boolean that indicates:
    • 1 = use line-drawing characters
    • 0 = use ASCII characters
  • intmat is an integer matrix representation of a maze
  • z is the drawing of the maze in ASCII or line-drawing characters

For example: pic←1 draw mat produces a line-drawing of the maze generated in the example above

 

      z←html intmat
  • intmat is an integer matrix representation of a maze
  • z is the HTML necessary to render the maze in a web browser. Save it to a native file and open the file in your browser.

For example: h←html mat produces an HTML representation of the maze generated in the example above

 

Solving the 2014 APL Problem Solving Competition – It’s All Right

This post is the continuation of the series where we examine some of the problems selected for the 2014 APL Problem Solving Competition.

The problems presented in Phase 1 of the competition were selected because they could be solved succinctly, generally in a single line of APL code. This makes them well suited for experimentation on TryAPL.org.

pythProblem 1 of Phase 1, entitled “It’s all right” reads,

“Write a dfn that takes the length of the legs of a triangle as its left argument, and the length of the hypotenuse as its right argument and returns 1 if the triangle is a right triangle, 0 otherwise.”

Test cases:
      3 4 {your_solution} 5
1
      2 3 {your_solution} 4
0

This uses the Pythagorean theorem – A² + B² = C². It’s trivial to implement an almost direct translation of this in APL – in a dfn, using ⍺[1] for A, ⍺[2] for B and for C yields:

right←{((⍺[1]*2)+(⍺[2]*2))=⍵*2}

This seems rather clunky though… what if we consider the problem as “Are the sums of the squares of each argument equal?” To get the sum of the squares, first we can use the composition *∘2 (composing the power function * with the constant 2) to mean “square” and +/ to mean “sum”, and combine them in a 2-item function train (also known as an “atop”): ((+/)*∘2)

then apply this to each argument:   ((+/)*∘2)¨⍺ ⍵

and compare the results for equality, resulting in the dfn:

right1←{=/((+/)*∘2)¨⍺ ⍵}

Still this seems clunky to me. Let’s see…squaring a number is just multiplying the number by itself, so, if we use the monadic commute operator with multiplication,   ×⍨
we get the equivalent of squaring. Then we can use that function in an inner product with addition to also get “sum of the squares”:   +.×⍨

The rest is essentially the same as in right1 above:

right2←{=/+.×⍨¨⍺ ⍵}

All three of these solutions solve the contest problem. The first one, right, is not commutative though – the legs of the triangle must be supplied as the left argument. However, right1 and right2 are commutative and the order of arguments is not significant.

Solving the 2014 APL Problem Solving Competition – it’s as easy as 1 1 2 3…

Competition LogoThe winners of the 2014 APL Problem Solving Competition were recognized at the Dyalog ’14 user meeting and had a chance to share their competition experiences. Emil Bremer Orloff won the student competition and received $2500 USD and an expenses-paid trip to Dyalog ’14, while Iryna Pashenkovska took first place among the non-student entries and received a complimentary registration to Dyalog ’14. Our congratulations to them for their fine efforts! This post is the first of a series where we will examine some of the problems selected for this year’s competition. The problem we’ll discuss first is Problem 3 from Phase I dealing with the Fibonacci sequence.

Problem 3 – Tell a Fib Write a dfn that takes an integer right argument and returns that number of terms in the Fibonacci sequence. Test cases:

      {your_solution} 10
 1 1 2 3 5 8 13 21 34 55
      {your_solution} 1
 1
      {your_solution} 0   ⍝ should return an empty vector

Essentially, you start with the sequence 1 1 and concatenate the sum of the last 2 numbers to the end and repeat until the sequence has the correct number of terms. In Python, a solution might look something like this:

def fib(n): # return the first n elements of the Fibonacci sequence
    result = []
    a, b = 0, 1
    while 0 < n:
        result.append(b)
        a, b = b, a+b
        n -= 1
    return result

You can write nearly the same code in APL:

r←fibLooping n;a;b
  r←⍬
  (a b)←0 1
:While 0<n
  r,←b
  (a b)←b,a+b
  n-←1
:EndWhile

While it’s possible to write APL code that looks like Python (or many other languages), one of the judging criteria for the competition is the effective use of APL syntax – in this case, by avoiding the explicit loop. Two ways to do this are 1) using recursion and 2) using the power operator ().

Recursive Solution

      fibRecursive←{⍵<3:⍵⍴1 ⋄ {⍵,+/¯2↑⍵}∇⍵-1}

The neat thing about recursion is that the function keeps calling itself with a “simpler” argument until it reaches a base case and then passes the results back up the stack. Here the base case occurs when the argument () is less than 3 and the function returns ⍵⍴1 – it’s either an empty vector, 1, or 1 1 for values of of 0, 1 and 2 respectively. It’s easy to illustrate the recursive process by adding some code (⎕←) to display information at each level of recursion.

      fibRecursive←{⍵<3:⎕←⍵⍴1 ⋄ {⎕←⍵,+/¯2↑⍵}∇⎕←⍵-1}
      fibRecursive 10
9
8
7
6
5
4
3
2
1 1
1 1 2
1 1 2 3
1 1 2 3 5
1 1 2 3 5 8
1 1 2 3 5 8 13
1 1 2 3 5 8 13 21
1 1 2 3 5 8 13 21 34
1 1 2 3 5 8 13 21 34 55

The recursive solution above is termed “stack recursive” in that it creates a new level on the function call stack for each recursive call. We can modify the code slightly to implement a “tail recursive” solution which Dyalog can detect and optimize by avoiding creating those additional function call stack levels. You can see the effect of this as each level computes its result immediately:

      fibTail←{⍺←0 1 ⋄ ⍵=0:⎕←1↓¯1↓⍺ ⋄ (⎕←⍺,+/¯2↑⍺)∇ ⎕←⍵-1}
      fibTail 10
        fibTail 10
9
0 1 1
8
0 1 1 2
7
0 1 1 2 3
6
0 1 1 2 3 5
5
0 1 1 2 3 5 8
4
0 1 1 2 3 5 8 13
3
0 1 1 2 3 5 8 13 21
2
0 1 1 2 3 5 8 13 21 34
1
0 1 1 2 3 5 8 13 21 34 55
0
0 1 1 2 3 5 8 13 21 34 55 89
1 1 2 3 5 8 13 21 34 55

Power Operator Solution

      fibPower←{⍵⍴({⍵,+/¯2↑⍵}⍣(0⌈⍵-1))1}

The power operator is defined as {R}←{X}(f⍣g)Y. When the right operand g is a numeric integer scalar – in this case (0⌈⍵-1), the power operator applies its left operand function f{⍵,+/¯2↑⍵} cumulatively g times to its argument Y, which in this case is 1. Similar to the recursive solution, we can add ⎕← to see what happens at each application:

      fibPower←{⍵⍴({⎕←⍵,+/¯2↑⍵}⍣(0⌈⍵-1))1}
      fibPower 10
1 1
1 1 2
1 1 2 3
1 1 2 3 5
1 1 2 3 5 8
1 1 2 3 5 8 13
1 1 2 3 5 8 13 21
1 1 2 3 5 8 13 21 34
1 1 2 3 5 8 13 21 34 55
1 1 2 3 5 8 13 21 34 55

In discussing this blog post with my fellow Dyalogers, Roger Hui showed me a rather neat construct:

      fibRoger←{1∧+∘÷\⍵⍴1}     ⍝ Roger's original version
      fibBrian←{1∧÷(+∘÷)\⍵⍴1}  ⍝ tweaked to make result match contest examples
      fibRoger 10
1 2 3 5 8 13 21 34 55 89
      fibBrian 10
1 1 2 3 5 8 13 21 34 55

I’ve added the parentheses around +∘÷ to make it a bit clearer what the left operand to the scan operator \ is. Let’s break the code down…

      ⍵⍴1 ⍝ produces a vector of ⍵ 1's 
1 1 1 1 1 ⍝ for ⍵=5

Then we apply scan with the function +∘÷, which in this case has the effect of adding 1 to the reciprocal of the previous result:

      (+∘÷)\1 1 1 1 1
1 2 1.5 1.666666667 1.6

Note that the denominator of each term is the corresponding element of the Fibonacci series…

1 ÷ 1 = 1
2 ÷ 1 = 2
3 ÷ 2 = 1.5
5 ÷ 3 = 1.6666666667
8 ÷ 5 = 1.6

To extract them, take the reciprocal and then the LCM (lowest common multiple) with 1.

      1∧÷1 2 1.5 1.666666667 1.6
1 1 2 3 5

What happens if you want to accurately compute larger Fibonacci numbers? There’s a limit to the precision based on the internal number representations. Because fibRoger and fibBrian delve into floating point representation, they’re the most limited. (Roger points out that the J language has native support for extended precision and does not suffer from this limitation.)

Dyalog has a system variable, ⎕FR, that sets the how floating-point operations are performed using either IEEE 754 64-bit floating-point operations or IEEE 754-2008 128-bit decimal floating-point operation. For most applications, 64-bit operations are perfectly acceptable, but in applications that demand a very high degree of precision, 128-bit operations can be used, albeit at some performance degradation. Using 64-bit operations, fibBrian loses accuracy after the 35th term while using 128-bits extends the accuracy to 68 terms.

Even setting ⎕FR to use 128-bit operations and the print precision, ⎕PP, to its maximum, we can only compute up to the 164th element 34-digit 8404037832974134882743767626780173 before being forced into scientific notation.

Can we go further easily? Yes we can! Using Dyalog for Microsoft Windows’ .NET integration, we can seamlessly make use of the BigInteger class in the System.Numerics .NET namespace.
First we need to specify the namespace and where to look for it…

      ⎕USING←'System.Numerics,system.numerics.dll'

Then we make a trivial change to our code to use BigInteger instead of native data types…

      fibRecursive←{⍵<3:⍵⍴⎕NEW BigInteger 1 ⋄ {⍵,+/¯2↑⍵}∇ ⍵-1}

And we’re done!

So the next time someone walks up to you and asks “What can you tell me about the 1,000th element of the Fibonacci series?”, you can confidently reply that it has 209 digits and a value of
43466557686937456435688527675040625802564660517371780402481729
08953655541794905189040387984007925516929592259308032263477520
96896232398733224711616429964409065331879382989696499285160037
04476137795166849228875.

For other APL Fibonacci implementations, check out this page.

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.