In a previous blog post I showed how a 2D convolver can be written in F# using Microsoft’s Accelerator system for GPGPU and SIMD x64 multicore programming targets. This blog article explains how the same 2D convolver can also be nicely expressed in C# and it illustrates a few cool C# and LINQ features including type inference, extension methods, the indexed Select operator and the Aggregate operator. I also give a more old school C# implementation at the end of this blog article. First, here is the code:

using System;
using System.Linq;

using Microsoft.ParallelArrays;

namespace AcceleratorSamples
{
    static partial class Convolver2D
    {
        static FloatParallelArray convolve(this FloatParallelArray a,
                                           Func<int, int[]> shifts, float[] kernel)
        {
            return kernel
                .Select((k, i) => k * ParallelArrays.Shift(a, shifts(i)))
                .Aggregate((a1, a2) => a1 + a2);
        }

        static FloatParallelArray convolveXY(this FloatParallelArray input, float[] kernel)
        {
            return input
                .convolve(i => new[] { -i, 0 }, kernel)
                .convolve(i => new[] { 0, -i }, kernel);
        }

        static void Main(string[] args)
        {
            const int inputSize = 10;

            var random = new Random(42);

            var inputData = new float[inputSize, inputSize];
            for (int row = 0; row < inputSize; row++)
                for (int col = 0; col < inputSize; col++)
                    inputData[row, col] = (float)random.NextDouble() * random.Next(1, 100);

            var testKernel = new[] { 2F, 5, 7, 4, 3 };

            var dx9Target = new DX9Target();
            var inputArray = new FloatParallelArray(inputData);
            var result = dx9Target.ToArray2D(inputArray.convolveXY(testKernel));

            for (int row = 0; row < inputSize; row++)
            {
                for (int col = 0; col < inputSize; col++)
                    Console.Write("{0} ", result[row, col]);
                Console.WriteLine();
            }
        }
    }
}

This version of the code was suggested by Mads Torgersen. I won’t explain the details of how 2D convolution is implemented here and I kindly ask you to first look over the original F# blog article. To recap a little, a 1D convolver computes a weighted average around each element of an array (or stream):

A 2D convolver works by first performing a convolution the x-direction and then performing a convolution in the y-direction. An x-direction convolution is shown below.

There are two basic operations that we need to express here:

• A multiplication of a shifted 2D input matrix by the i-th element of a kernel 1D matrix. This is accomplished with an indexed Select operator.
• The addition of all the elements in a row (or column) of a 2D matrix. This is accomplished with an Aggregate operator.

First let me explain what C# extension methods are by using the code below as a demonstration.

using System;

namespace Extension_Methods_Demonstration
{
    public static class Extension_Methods_Experiment
    {

        public static int DoubleFunction(int i)
        {
            return 2 * i;
        }

        public static int Double(this int i)
        {
            return 2 * i;
        }

        static void Main(string[] args)
        {
            int x = 12;
            Console.WriteLine(DoubleFunction(x));
            Console.WriteLine(x.Double());
        }
    }
}

Here we have two ways of encapsulating a computation that doubles an int value. The first approach uses a regular static method method DoubleFunction which takes an int value as an argument and returns an int value as a result i.e. something that behaves like a function. It might be nice to also express the doubling operation as a regular method invocation on variables of type int. However, we don’t always have access to the definition for all the types we want to extend in this way. Extension methods allow us to express functions over types even when we don’t have access to the definition of the type. The syntax is very similar to the syntax used for static methods however the type of the object being extended is indicated as the first parameter and is identified with the this keyword. The program above illustrates both the static method style invocation of a double operation as well as the invocation of a Double method on an variable of type int.

Before we explain what an indexed Select operator does we first demonstrate what a regular Select operator does:

using System;
using System.Linq;

class Select_Demonstration1
{
    static void Main(string[] args)
    {
        int[] numbers = { 5, 4, 1, 3, 9, 8, 6, 7, 2, 0 };
        var y = numbers.Select(n => 2 * n);
        foreach (var i in y)
            Console.Write(i + " ");
        Console.WriteLine();
    }
}

When run this program produces the output:

10 8 2 6 18 16 12 14 4 0

Looking at the example inside out, let us focus on the lambda expression n => 2 * n first. This defines an anonymous function (static method) which takes an int argument n and returns an int result which is twice the value of the input value. It is short-hand for:

delegate (int n) { return 2 * n; }

except the lambda expression is a little shorter and all the types are inferred. The lambda expression form is inspired by lambda expressions in the lambda calculus.

The code snippet above uses the Select operator to apply the same operation (function) to event element of an array to yield a new collection. The operation to be performed is specified here using a lambda expression which just doubles its input. This is very much like the higher order map function in functional languages e.g. List.map in F#. The result of applying the Select operator is a variable y that implements the IEnumerable interface which allows us to use this variable in foreach loops. Note that we have used the var keyword to indicate that we wish the type of y to be inferred. If you are a functional programmer you might assume that Select is exactly like map then you would expect the type of y to be int[] but sadly you would be wrong (but let’s not get bogged down with that here). The overloaded type of Select used here is:

public static IEnumerable<TResult> Select<TSource, TResult>  (this IEnumerable<TSource> source,    Func<TSource, TResult> selector);

Next we turn out attention to an extension method called indexed Select. This operator can be used to apply some operation to every element of a collection (e.g. an array) but unlike the basic Select this operator also provides the index value of each element. The argument to select is a function that takes a pair that contains an element and its index taken from the source array. The result of the function specifies the value at the corresponding index in the output collection. Here’s an example that multiplies each value in an array by its index value:

using System;
using System.Linq;

class Select_Demonstration2
{
    static void Main(string[] args)
    {
        int[] numbers = { 5, 4, 1, 3, 9, 8, 6, 7, 2, 0 };
        var y = numbers.Select((n, i) => n * i);
        foreach (var i in y)
            Console.Write(i + " ");
        Console.WriteLine();
    }
}

which produces the output:

0 4 2 9 36 40 36 49 16 0

i.e. 5*0, 4*1, 1*2, 3*3, 9*4 etc.

This version of Select has the overloaded type:

static IEnumerable<TResult> Select<TSource, TResult>  (this IEnumerable<TSource> source,    Func<TSource, int, TResult> selector);

Now we can begin to understand what the Select portion of the convolve method does:

static FloatParallelArray convolve(this FloatParallelArray a,
                                   Func<int, int[]> shifts, float[] kernel)
{
    return kernel
        .Select((k, i) => k * ParallelArrays.Shift(a, shifts(i)))
        .Aggregate((a1, a2) => a1 + a2);
}

First we see that convolve is an extension method because the first parameter is preceded by the this keyword. In particular this method extends the FloatParallelArray type with this convolve method. The convolve method takes two arguments:

• Func<int, int[]> shifts: this function takes as input an int value which specified the index value of a kernel element and it returns a two element array which specifies the shifts in the x-direction and the y-direction.
• float[] kernel: this is the 1D kernel array

We now see that the Select operator multiplies each kernel value by the input array shifted according to some function that is specified as an argument to convolve. Let’s consider a specific invocation:

kernel.Select((k, i) => k * ParallelArrays.Shift(a, i => new[] { 0, –i }));

Here we use a lambda expression to indicate the shift degree in the x-direction should match the index value of the current kernel element and the shift degree in the y-direction is always zero. This helps to implement a portion of a 1D convolution of a 2D input array as shown below.

To complete the 1D convolution of the resulting collection of 2D arrays we need to add these arrays together in a point-wise manner. We can achieve this with the extension method Aggregate:

public static TSource Aggregate<TSource>  (this IEnumerable<TSource> source, Func<TSource, TSource, TSource> func);

The Aggregate extension method is very similar to a left fold in functional languages e.g. List.fold in F#. Informally, this method takes a collection of some values e.g. {A, B, C, D} and combines them with some operator, say ͦ, by computing ((A ͦB) ͦC) ͦ D). For example, we can specify the + operation to add up a collection of numbers {A, B, C, D} using Aggregate  which computes A + B + C + D. So to compute the sum of the collection of 2D arrays produced by the Selectoperator we simply need to apply the opertator with a lambda expression that adds two FloatParallelArrays. We simply used the + operator which is overloaded for FloatParallelArray.An example of the application of Aggregate to the result of the collection produced by the previous Select is shown below.

This completes the description of the convolve extension method which implements a 1D convolution of a 2D input array. Now we define a 2D convolution convolveXY extension method by specifying a convolution in the y-direction followed by a convolution in the x-direction:

static FloatParallelArray convolveXY(this FloatParallelArray input, float[] kernel)
{
    return input
        .convolve(i => new[] { -i, 0 }, kernel)
        .convolve(i => new[] { 0, -i }, kernel);
}

Note the use of anonymous array initializers i.e. the type of the two element arrays is not explicitly specified.

The Main method uses fairly straight forward C# code to create some test input data, create a DX9 Accelerator target and run and display the results. The code makes extensive use of var to infer the types of variables.

One advantage of making convolve and convolveXY into extension methods it that we can use the nice fluent style of composition which allows us to chain together transformations over collections of FloatParallelArrays.

For comparison, here is another more straight forward C# implementation of the convolver.

using System

using Microsoft.ParallelArrays;

namespace AcceleratorSamples
{
    partial class Convolver2D
    {
        static FloatParallelArray convolve(Func<int, int[]> shifts, float[] kernel,                                            int i, FloatParallelArray a)
        {
            FloatParallelArray e = kernel[i] * ParallelArrays.Shift(a, shifts(i));
            if (i == 0)
                return e;
            else
                return e + convolve(shifts, kernel, i - 1, a);
        }

        static FloatParallelArray convolveXY(float[] kernel, FloatParallelArray input)
        {
            FloatParallelArray convolveX                  = convolve(i => new int[] { -i, 0 }, kernel, kernel.Length - 1, input);
            return convolve(i => new int[] { 0, -i }, kernel, kernel.Length - 1, convolveX);
        }

        static void Main(string[] args)
        {

            int inputSize = 10;

            Random random = new Random(42);

            float[,] inputData = new float[inputSize, inputSize];
            for (int row = 0; row < inputSize; row++)
                for (int col = 0; col < inputSize; col++)
                    inputData[row, col] = (float)random.NextDouble() * random.Next(1, 100);

            float[] testKernel = new float[]{2, 5, 7, 4, 3} ;

            DX9Target dx9Target = new DX9Target();
            FloatParallelArray inputArray = new FloatParallelArray(inputData);
            float[,] result = dx9Target.ToArray2D(convolveXY (testKernel, inputArray));

            for (int row = 0; row < inputSize; row++)
            {
                for (int col = 0; col < inputSize; col++)
                    Console.Write("{0} ", result[row, col]);
                Console.WriteLine();
            }
        }
    }
}

This blog explains how Escher‘s fish wood carving can be expressed in F# using functional geometry which was first presented by Peter Henderson in 1982. This very functional and compositional style of describing pictures is nicely expressed in F#. Functional geometry provides us with a nice algebra for pictures which has been used by the author as the basis of a system for automatically computing the layout of digital circuits expressed in functional languages (e.g. Lava).

First, here is the picture:

And here are the core F# definitions that are used to make the drawing:

The rest of this blog explains how these functions are used to draw the Escher fish picture. A distinctive aspect of drawing pictures using functional geometry is compositional nature in which a composite picture is computed from component pictures by a series of transformations like rotations and scaling operations. In particular there is a lack of concrete coordinates because pictures are typically composed using higher order combinators that work with relative locations e.g. “picture A is to the left of picture B” or “picture A is below picture B”.

This picture is produce by transforming four basic tiles called p, q, r and s:

p                                     q                                    r                                    s

The functions used to compute the compute the final drawing operate over a record type that describes points in two dimensions along with some associated members for certain kinds of point operations

Note that the ~ symbol is used to specify the overloading of a unary operation (negation in this case). We could define many other overloaded members for *, + etc. but here we define just the operations we need for the functional geometry system. We also use the Point type to describe two-element vectors.

We represent a picture as a function which takes three vectors a, b, and c and produces a drawing with its origin translated by the vector a and with its horizontal and vertical dimensions bounded by the vectors b and c as illustrated below.

A picture is function of the values of the a, b and c vectors which are then used to produce a drawing as a side effect.  A picture in our functional geometry system has the type:

To draw a line we need a graphics canvas and a pen plus a function plotLine which takes the following parameters:

·         m and n which specify the size of the picture grid

·         a, b, c which are the vectors which specify the origin and the orientation

·         p1 and p2 which specify two Point values which identify the stat and end points of a line segment.

A picture can be expressed as a collection of line segments by using the grid function:

The grid function is defined as follows:

This function works by simply applying the plotLine function to each line segment with the appropriate grid size, origin and orientation. The ignore function is used to ensure the return type of this function is unit.

These basic p, q, r and s tiles have little regularity so we define the contents of these tiles by giving their bounding boxes and a list of line segments. Here is the definition of tile p:

The pointfy function is used to convert the list of pairs of integers to pairs containing Points.

A useful picture is the empty picture which contains no line segments so we anoint it with a special name:

This picture generating function is more useful than you might think as we shall shortly demonstrate.

A useful higher order function is rot which rotates a picture by 90 degrees.

This is a higher order function because it takes a picture (a function) as its input and returns a picture (a function) as its result. We call such picture to picture functions picture combinators. The application of rot to the tile p is demonstrated below.

                  p                               rot p

We will need to place one picture beside another so we define a combinator which performs this task:

Note that this definition makes use of the overloaded members +, * and – for the type Point. The picture of p beside q is shown below.

Similarly we need to place one picture above another picture:

We define t to be a quartet formed by placing the four basic riles in a square.

Here is the picture of t which shows how nicely the four basic fish tiles fit together:

Another useful picture combinator is cycle which creates a new picture by making successive rotations of an input picture and laying them out using the quartet pattern.

A picture of cycle (rot t) is shown below.

We now define two pictures which have a very interesting property.

                             side1                                                 side2

The amazing property of the picture t is that when reduced to half its size, it sits perfectly on top of itself (in two different ways)!

We are close to making the final fish picture but first we must define a function that can take nine pictures as its input and produce a regular three by three arrangement.

This picture combinator can be used to define a picture called corner.

which is shown below.

The final fish picture, called Square Limit, can now be defined by applying cycle to this corner.

This produces the picture shown at the top of this blog.

The original deconstruction of Escher’s fish was reported by J. L. Locher in 1971 in the book “The World of M. C. Escher”. Here we present a “reconstruction” using a nice algebra of picture combinators. For more details about how the fish are drawing using functional geometry look at Henderson’s papers e.g. Functional Geometry, or Functional Geometry, or Functional Geometry (yes, they really are all called the same but are somewhat different). The original illustrations for the 1982 paper were done by Mary Sheeran using a functional geometry package implemented in UCSD Pascal. Note that the origin for the drawings on this blog is at the top right hand corner although the origin for the drawings in the original paper is at the bottom right.

Many others have re-implemented Henderson’s fish. Here are just a few of them:

·         Frank Buss http://www.frank-buss.de/lisp/functional.html

·         Bill Clementson http://bc.tech.coop/blog/050213.html

PanBrowser incorporates Conal Elliott‘s ideas on functional images into an F# library.

I believe the idea behind functional geometry can be applied to description of many sophisticated entities including pictures, circuits and music. In each situation we have to try and understand the underlying structures and look for regular patterns and hierarchy and powerful ways of combining sub-creations into composite creations.

The complete F# program for drawing Escher’s fish is shown below which I have tested with F# in Visual Studio 2010 Beta 2.

Satnam Singh
http://research.microsoft.com/~satnams