A small programming exercise in APL #2: Combinations

In this post I'm going to continue my APL exploration by working in a possible solution for the problem of generating all the combinations of 'n' elements of an array.

Generating the 'n' combinations of an array means generating a sequence of all possible 'n' unique elements from the original array . For example, given the array 12, 34, 35, 65 all possible '2' combinations of this array are:

• 12, 34
• 12, 35
• 12, 65
• 34, 35
• 34, 65
• 35, 65

Notice that order is not important. That is "12, 34" is considered to be the same as "34, 12". Generating all 'n' permutations of the elements of an array may be an interesting problem for a future post.

One of my goals was to be as idiomatic as possible (with my limited APL knowledge!). Because of this, I will try avoid using explicit loops or conditionals and instead use array operators.

Strategy

The general strategy for solving this problem was to calculate all possible boolean arrays having 'n' bits and use Compress to extract the elements.

For example, in the array 12, 34, 35, 65 all possible boolean vectors having two bits on are:

• 1 1 0 0
• 1 0 1 0
• 1 0 0 1
• 0 1 1 0
• 0 1 0 1
• 0 0 1 1

Using these vectors in APL we get the elements we need:

      a
12 34 35 65
      1 1 0 0 / a
12 34
      1 0 1 0 / a
12 35
      1 0 0 1 / a
12 65
      0 1 1 0 / a
34 35
      0 1 0 1 / a
34 65
      0 0 1 1 / a
35 65

This strategy is not efficient in space or time but it allowed me to explore a solution without using imperative constructs or recursion.

Generating the boolean arrays

To generate the boolean arrays we start by generating all integer values required to encode n-bit values:

      a ← 12 34 35 65
      ⍳ ( ¯1 + 2 * ⍴ a)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

When converting these numbers to binary we get all possible boolean vectors for 4-bit values. This is same as the same as the length of our sample array. We can see these arrays if we use Encode and Each:

     { ((⍴ a) ⍴ 2) ⊤ ⍵ } ¨ ⍳ ( ¯1 + 2 * ⍴ a)
 0 0 0 1  0 0 1 0  0 0 1 1  0 1 0 0  0 1 0 1  0 1 1 0  0 1 1 1  1 0 0 0\
  1 0 0 1  1 0 1 0  1 0 1 1  1 1 0 0  1 1 0 1  1 1 1 0  1 1 1 1

I can reshape this array to get a better representation:

      15 1 ⍴ { ((⍴ a) ⍴ 2) ⊤ ⍵ } ¨ ⍳ ( ¯1 + 2 * ⍴ a)
 0 0 0 1
 0 0 1 0
 0 0 1 1
 0 1 0 0
 0 1 0 1
 0 1 1 0
 0 1 1 1
 1 0 0 0
 1 0 0 1
 1 0 1 0
 1 0 1 1
 1 1 0 0
 1 1 0 1
 1 1 1 0
 1 1 1 1

Now we are only interested in boolean arrays with 'n' number of '1' bits. For n = 2 we can get only those elements using the Compress operator:

      seq ← ⍳ ( ¯1 + 2 * ⍴ a)
      n ← 2
      ({ n =  +/ ((⍴ a) ⍴ 2) ⊤ ⍵} ¨ seq) / seq
3 5 6 9 10 12

Again we can visualize these values using reshape and encode:

      6 1 ⍴ { ((⍴ a) ⍴ 2) ⊤ ⍵ } ¨ ({ n =  +/ ((⍴ a) ⍴ 2) ⊤ ⍵} ¨ seq) / seq
 0 0 1 1
 0 1 0 1
 0 1 1 0
 1 0 0 1
 1 0 1 0
 1 1 0 0
 

Selecting desired elements

Now that we have our boolean arrays we can pick all the values using Compress:

     {(((⍴ a) ⍴ 2) ⊤ ⍵)/a} ¨ ({ n =  +/ ((⍴ a) ⍴ 2) ⊤ ⍵} ¨ seq) / seq
 35 65  34 65  34 35  12 65  12 35  12 34

And again we can reshape this array to see it better:

      6 1 ⍴ {(((⍴ a) ⍴ 2) ⊤ ⍵)/a} ¨ ({ n =  +/ ((⍴ a) ⍴ 2) ⊤ ⍵} ¨ seq) / seq
 35 65
 34 65
 34 35
 12 65
 12 35
 12 34

Finally we can pack these expressions in a function:

∇r←a ncombinations1 n;seq
  seq ← ⍳ ( ¯1 + 2 * ⍴ a)
  r ← {(((⍴ a) ⍴ 2) ⊤ ⍵)/a} ¨ (({ n =  +/ ((⍴ a) ⍴ 2) ⊤ ⍵} ¨ seq) / seq\
)
∇

We can use it like this:

      1 2 3 4 5 ncombinations1 3
 3 4 5  2 4 5  2 3 5  2 3 4  1 4 5  1 3 5  1 3 4  1 2 5  1 2 4  1 2 3

APL has a build-in Binomal operator which allows us to calculate the number of combinations. For example:

      a ← 1 2 3 4 5
      (3 ! ⍴ a) = ⍴ a ncombinations3 3
1

Final words

It was very interesting (and difficult) to try to write a solution of this problem using only array operations. Reading the definion of the ncombinations1 I noticed that there are many expressions nested in parenthesis (maybe related to my limited APL knowledge). I know that APL developers value terseness, so I tried to reduce the size of the expressions. Here is a new version:

∇r←a ncombinations3 n;seq
  s ← ⍳ ( ¯1 + 2 * ⍴ a)
  b ← 2⍴⍨⍴a
  r ← {a/⍨b⊤⍵}¨s/⍨{n=+/b⊤⍵}¨s
∇

I was able to remove a couple of parenthesis by taking advantage of the Commute operator.

GNU APL was used to create the examples in this post.

A small programming exercise in APL #1

This post shows a possible solution written in APL for the following programming problem:

Given an array, determine if a value 'number' is found in every consecutive segment of 'size' elements.

For example, given this array:

1,2,1,3,4,1

This predicate is true for number = 1 and size = 2. Because '1' is found in (1,2), (1,3) and (4,1).

A possible APL solution

This is a possible solution to this problem (using my limited APL knowledge).

∇ z←array arraysegments args
  number ← args[1]
  size ← args[2]
  z ← ∧/ ({number∊⍵} ⍤1) (((⍴ array) ÷ size) , size) ⍴ array
∇

This function is called arraysegments and is a dyadic function. It receives the array to validate as the left argument and on the right argument a two-value array with the value ('number') to find inside the array segment and the length of the segments.

An example of using this function

      1 2 1 3 4 1 arraysegments 1 2
1
      1 2 1 3 4 1 arraysegments 2 2
0
      1 2 1 3 4 1 arraysegments 1 3
1
      300 200 300 300 400 500 arraysegments 300 3
1
      300 200 300 300 400 500 arraysegments 400 3
0

To see how this function works I'm going to start by defining some sample data:

      array ← 1 2 1 3 4 1 
      number ← 1
      size ← 3

1. First, start by reshaping the input array into a matrix with each row being

      (((⍴ array) ÷ size) , size)
2 3      
      (((⍴ array) ÷ size) , size) ⍴ array
1 2 1
3 4 1

2. Now that we have the 'groups' as the rows of the new matrix, we apply the membership operation to each element of the matrix. To do this, we use the rank operator in conjunction with an inline function to test the membership:

      ({number∊⍵} ⍤1) (((⍴ array) ÷ size) , size) ⍴ array
1 1

As a final step we apply the reduce operator '/' with the And operator to see if the input number exists in every group:

      ∧/ ({number∊⍵} ⍤1) (((⍴ array) ÷ size) , size) ⍴ array
1

There are some problem with this solution. For example if the array cannot be splitted in groups of the specified size you get a runtime error.

Reading APL #1: Roman numerals to decimal

As part of a new recurring series I’m going to take a small APL function or program and try to break it down in to its parts. I hope this is going to help me get a better understanding of the language.

Reading APL

This article: The APL Programming Language Source Code has a nice introduction to the language and its history. One sentence in this article is key to understand how to read APL code:

Order of evaluation: Expressions in APL are evaluated right-to-left, and there is no hierarchy of function precedence

Two additional concepts are useful to understand when trying to read code:

1. A monadic function : an operation applied to the right expression
2. A dyadic function : a function applied to two arguments ( left and right )

One example of these elements is the following:

      2 2 ⍴⍳⍴ 8 8 8 8
1 2
3 4

Here is an example of evaluating this expression from right to left.

      myarray←8 8 8 8
      myarray
8 8 8 8
      length_of_my_array← ⍴ myarray
      length_of_my_array
4
      new_array← ⍳ length_of_my_array
      new_array
1 2 3 4
      2 2 ⍴ new_array
1 2
3 4

Note that here we use ⍴ both as a monadic function (Shape) and as a dyadic function (Reshape).

The example

For this post I’m going to use a small function I found in the APL Utils repository by Blake McBride . This function takes a string representing a roman number and converts it back to an integer:

∇z←Romanu a
 z←+/(¯1+2×z≥1↓z,0)×z←(1 5 10 50 100 500 1000)['IVXLCDM'⍳a]
∇

An example of executing this function:

      Romanu 'XVII'
17
      Romanu 'IX'
9

Reading the code

Let’s start by reading the function definition:

∇z←Romanu a
   ...
∇

This code represents the definition of the Romanu function with an a argument.

Now, the interesting part is reading the body of the function which performs the actual calculations.

z←+/(¯1+2×z≥1↓z,0)×z←(1 5 10 50 100 500 1000)['IVXLCDM'⍳a]

As described in the previous section we need to read this code from right to left. To do this we need to ‘parse’ it and extract the largest complete expression from the right:

(1 5 10 50 100 500 1000)['IVXLCDM'⍳a]

This expression is a particular example of an bracked indexing expression(https://aplwiki.com/wiki/Bracket_indexing). A simpler example of this is:

      (10 20 30)[1]
10
      (10 20 30)[2]
20
      (10 20 30)[3]
30
      (10 20 30)[3 1 1 2]
30 10 10 20
      a←10 20 30
      a[3 1 1 2]
30 10 10 20

As shown above this expression is used to select elements of an array into a new array.

In the Romanu function this expression is used to assign a value to eacho of the roman numeral 'digits':

'IVXLCDM'⍳a

Here the dyadic form of ⍳ or ‘Index of’(https://aplwiki.com/wiki/Index_Of) . This expression returns the indices of the left side array of the elements on the right. Here are some examples of this expression in action:

      (10 20 30) ⍳ (20 20 10 20 30 10)
2 2 1 2 3 1

When using this expression with the roman digits string we get the equivalent positions:

      'IVXLCDM' ⍳ 'XVII'
3 2 1 1

As shown above we get the positions in the roman digits string. Returning to our original expression we use the selection expression to get the value of each roman digit in decimal form:

      (1 5 10 50 100 500 1000)['IVXLCDM' ⍳ 'XVII']
10 5 1 1

As you can see here we almost have the result for the conversion of XVII. We just have to sum these numbers up to get 17.

Going back to the example and continuing reading for right to left we find and assignment to the z variable.

z←(1 5 10 50 100 500 1000)['IVXLCDM'⍳a]

The rest of the expression is going to take care of subtracting the value required for converting numbers like: ‘IX’.

Reading the expression from right to left we find:

(¯1+2×z≥1↓z,0)×z←(1 5 10 50 100 500 1000)['IVXLCDM'⍳a]

This expression is a multiplication (×). For arrays it applies to operation element wise :

      2 3 × 4 5
8 15

The interesting part here is the expression being used for the left operand of the multiplication:

(¯1+2×z≥1↓z,0)

Again we are going to analyze this expression from right to left.

z,0

This expression returns the an array with a zero at the end:

     (10 20 30),0
10 20 30 0
      z←(1 5 10 50 100 500 1000)['IVXLCDM' ⍳ 'IX']
      z,0
1 10 0

Then the following operation drops the first element of the array:

      1↓z,0
10 0

The drop expression returns an array without the first ‘n’ elements:

      3 ↓ 20 30 40 30 20 10
30 20 10

The following section is very interesting . We use the greater equal expression ≥ to determine which elements of the array are greater or equal than the second array:

       10 20 30 ≥ 5 21 30
1 0 1

Applying this to our example for ‘IX’:

      z←(1 5 10 50 100 500 1000)['IVXLCDM' ⍳ 'IX']
      z≥1↓z,0
0 1

This is interesting: we get an array that indicates which entries are lower than its successor. In the case of IX we are getting 0 1 saying that the first digit needs to be subtracted from the second. The final part performs this operation.

      z←(1 5 10 50 100 500 1000)['IVXLCDM' ⍳ 'IX']
      z≥1↓z,0
0 1
      ¯1+2×z≥1↓z,0
¯1 1
      (¯1+2×z≥1↓z,0)×z
¯1 10

Now we can finally sum all the numbers on the resulting array and get the decimal number. This is performed with the reduce operation / (https://aplwiki.com/wiki/Reduce) with the + operator.

      +/ 10 20 30
60
      10 + 20 + 30
60

Applying this operation to our partial expression returns the expected value:

      (¯1+2×z≥1↓z,0)×z
¯1 10

      +/(¯1+2×z≥1↓z,0)×z
9

We can see this process in action by executing all the parts with an interesting number like MMCXXIX (2129).

      z←(1 5 10 50 100 500 1000)['IVXLCDM'⍳a]
      z
1000 1000 100 10 10 1 10
      z,0
1000 1000 100 10 10 1 10 0
      1↓z,0
1000 100 10 10 1 10 0
      z≥1↓z,0
1 1 1 1 1 0 1
      ¯1+2×z≥1↓z,0
1 1 1 1 1 ¯1 1
      +/(¯1+2×z≥1↓z,0)×z
2129

A small experiment with fragment shaders

I wanted to work on an experiment that allowed me to learn a little bit about modern graphics programming with OpenGL (with shaders that is). A nice choice is the 'hello world; of graphics programming: rendering the Mandelbrot set.

In this post I'm going to record my experiences from zero to a basic rendering of the Mandelbrot set. Here is how the final product looks like:

1.1 The experiment

To get the "native" feel I decided to create the example as a C++ Windows program (but using mostly C). There are many options to do this, but the following combination worked for me:

• A C++ Compiler : Visual Studio build tools
• An IDE: VSCode C++ integration (https://code.visualstudio.com/docs/cpp/config-msvc)
• OpenGL libraries and loader: Windows SDK (https://developer.microsoft.com/en-us/windows/downloads/windows-10-sdk/) + GLAD (https://glad.dav1d.de/)
• A build system (although the experiment is quite small): CMake
• A simple GUI library: GFLW (https://www.glfw.org/)

I'm quite impressed with the VSCode C++ support .

1.2 The program

The structure of the program it's very simple since it is a hello world program. We start with some initialization boilerplate (quite short since GLFW abstracts away many windowing system details):

GLFWwindow *window;

glfwSetErrorCallback(error_callback);

if (!glfwInit())
{
  return -1;
}

window = glfwCreateWindow(
    800,
    600,
    "Shader example",
    NULL,
    NULL);

if (!window)
{
  glfwTerminate();
  return -1;
}
glfwMakeContextCurrent(window);
gladLoadGL();

This code creates a 800x600 window ready to use with OpenGL. To access OpenGL I'm using a loaded called GLAD.

1.3 Strategy for rendering the Mandelbrot set

In this experiment I am going to use the "escape time algorithm" to render the image of the Mandelbrot set. This algorithm is very simple and easy to implement using fragment shaders.

Before we start implementing this algorithm we are going to define geometry to and apply a fragment shader to it. Since we are going to render a 2D image, the easiest way to do this is to create two triangles that fill the viewport. Since OpenGL uses normalized coordinates we can define this triangles using values between -1.0 and 1.0 using the following code:

GLfloat vertices[] = {
    -1.0f, 1.0f,
    1.0f, 1.0f,
    1.0f, -1.0f,

    1.0f, -1.0f,
    -1.0f, -1.0f,
    -1.0f, 1.0f};

GLuint vertex_buffer;
glGenBuffers(1, &vertex_buffer);
glBindBuffer(GL_ARRAY_BUFFER, vertex_buffer);
glBufferData(
    GL_ARRAY_BUFFER,
    sizeof(vertices),
    vertices,
    GL_STATIC_DRAW);

This code is going to define two triangles that are highlighted here:

This code is also making these coordinates as the active ones (glBindBuffer).

Now we need both a vertex shader to apply transformations to the vertex data and a fragment shader to provide color to our pixels.

Our vertex shader is really simple since we are not performing transformations to the triangles:

#version 110
attribute vec2 pos;
void main() {
  gl_Position = vec4(pos, 0.0, 1.0); 
}

In this vertex shader we can manipulate the vertices of the triangles that we are going to draw. The X and Y values of the vertices are passed using the pos attribute. This is not automatic, we need to specify the data that is passed to the vertex shader in the C++ program for the pos attribute. This is accomplished by using the following code:

GLuint vpos_location;
vpos_location = glGetAttribLocation(program, "pos");
glEnableVertexAttribArray(vpos_location);
glVertexAttribPointer(
    vpos_location,
    2,
    GL_FLOAT,
    GL_FALSE,
    sizeof(float) * 2, 0);

One of the most interesting aspects of this snippet is the [glVertexAttribPointer](https://www.khronos.org/registry/OpenGL-Refpages/gl4/html/glVertexAttribPointer.xhtml). This function specifies the way the values are going to be extracted from the array data. Here we specify:

1. vpos_location the attribute that we are configuring
2. 2 the number of components (remember that pos is vec2)
3. GL_FLOAT the data type of the data element
4. GL_FALSE the data is not normalized. This seems to be important for integer data (here we use floating point data, more info here: https://gamedev.stackexchange.com/questions/10859/glvertexattribpointer-normalization).
5. sizeof(float) * 2 The offset between consecutive elements. This is useful when the array has different kinds of data elements. For example vertices and normals mixed in the same array. This parameter must be used to skip undesired data elements. In our case 0 is also a valid value since we only have vertex data.

The boilerplate code to compile this shader is the following:

GLint shader_compiled;
GLuint vertex_shader;
vertex_shader = glCreateShader(GL_VERTEX_SHADER);

glShaderSource(vertex_shader, 1, &vertex_shader_src, NULL);
glCompileShader(vertex_shader);

glGetShaderiv(vertex_shader, GL_COMPILE_STATUS, &shader_compiled);
if (shader_compiled != GL_TRUE)
{
  std::cout << "Vertex shader not compiled";
  GLchar message[1023];
  GLsizei log_size;
  glad_glGetShaderInfoLog(vertex_shader, 1024, &log_size, message);
  std::cout << message;
}

When you are starting with shaders, the call to glGetShaderiv is useful . This code returns error information that occurred when compiling the shader.

There is an important part of the process that is always confusing to me. Many things in the OpenGL API change a global state aspect of the functionality. For example the glBindBuffer function used above is going to determine the set of vertices used by the glDrawArrays function.

Here's the code of the fragment shader:

#version 110
uniform vec4 boundaries;  
void main()  {
   float x0, y0,x ,y;
   x0 = gl_FragCoord.x*((boundaries.z - boundaries.x)/800.0) + boundaries.x ;
   y0 = gl_FragCoord.y*((boundaries.w - boundaries.y)/600.0) + boundaries.y ;
   int maxIteration, iteration;
   maxIteration = 256;
   iteration = 0;
   while((x*x + y*y <= 4.0) && (iteration < maxIteration)) {
      float tmp;
      tmp = x*x - y*y + x0;
      y = 2.0*x*y + y0;
      x = tmp;
   iteration = iteration + 1;
   }
   gl_FragColor = vec4(vec3(0.0, float(iteration)/256.0, 0.0 ),1.0);
}

This code contains a simple implementation of the escape time algorithm described in Wikipedia here: https://en.wikipedia.org/wiki/Plotting_algorithms_for_the_Mandelbrot_set#Unoptimized_naïve_escape_time_algorithm .

For simplicity this shader chooses only between shades of green for the colors.

1.4 Zooming in

As part of this experiment I wanted to have the possibility to zoom a particular area of the Mandelbrot set. To add support for this we need to pass the top-left and bottom right coordinates of the area when want to render. This is accomplished by declaring a uniform to pass this value from the C++ program.

Here is the code that uses this parameter in the fragment shader:

uniform vec4 boundaries;  
...
   x0 = gl_FragCoord.x*((boundaries.z - boundaries.x)/800.0) + boundaries.x ;
   y0 = gl_FragCoord.y*((boundaries.w - boundaries.y)/600.0) + boundaries.y ;

And here is the code that specify the value of the boundaries uniform in the C++ program:

int coordinatesUniformLocation = glGetUniformLocation(program, "boundaries");
glUniform4f(coordinatesUniformLocation, boundaries.x1, boundaries.y1, boundaries.x2, boundaries.y2);

We can manipulate the values of the boundaries array using a mouse click handler like this:

glfwSetMouseButtonCallback(window, mouseCallback);
...
void mouseCallback(GLFWwindow* window, int button, int action, int mods) {
  if (action == GLFW_PRESS) {
     double xpos, ypos;
     glfwGetCursorPos(window, &xpos, &ypos);

ypos = 600 - ypos;
     double xclick, yclick;
     xclick = xpos*((boundaries.x2 - boundaries.x1)/800.0) + boundaries.x1;
     yclick = ypos*((boundaries.y2 - boundaries.y1)/600.0) + boundaries.y1;

     double currentWidth = (boundaries.x2 - boundaries.x1) - (boundaries.x2 - boundaries.x1)/10;
     double currentHeight = (boundaries.y2 - boundaries.y1) - (boundaries.y2 - boundaries.y1)/10;
     boundaries.x1 = (float)( xclick - currentWidth/2);
     boundaries.x2 = (float)( xclick + currentWidth/2);

     boundaries.y1 = (float)( yclick - currentHeight/2);
     boundaries.y2 = (float)( yclick + currentHeight/2);
   }
}

Finally we define the code draw loop like this:

while (!glfwWindowShouldClose(window))
{
  int width, height;
  glfwGetFramebufferSize(window, &width, &height);
  glViewport(0, 0, width, height);
  glClear(GL_COLOR_BUFFER_BIT);
  glUseProgram(program);
  int coordinatesUniformLocation = glGetUniformLocation(program, "boundaries");

  glUniform4f(coordinatesUniformLocation, boundaries.x1, boundaries.y1, boundaries.x2, boundaries.y2);

  glDrawArrays(GL_TRIANGLES, 0, 6);
  glfwSwapBuffers(window);
  glfwPollEvents();
}

1.5 Conclusion

This was a fun experiment!. I got a small glimpse on how to work with OpenGL. Future posts are going to explore even further and maybe use other programming languages to explore OpenGL.

Code for this experiment can be found here: https://github.com/ldfallas/openglmandelbrot .