Friday, February 25, 2011

Merge sequences in parallel

In practice you may find yourself in a situation when you have several sequences of data that you want to drain in parallel and merge the results into a single sequence. This is what Merge combinator for sequences does. Reactive Extensions had it for enumerable and observable sequences however at some point it was decided to no longer support it for enumerable sequences in Reactive Extensions as it may be expressed through the same combinator for observable sequences. 

I find it quite interesting to implement Merge combinator for enumerable sequences but to make it more interesting let’s change behavior a little bit. The original behavior of the EnumerableEx.Merge was to drain source sequences as fast as workers can and buffer results in a resulting queue until it is consumed by merged sequence enumerator. It can be implemented (maybe this is not the most efficient way but still) using

  • BlockingCollection<T> to synchronize producers that drain original sequences and consumer that merges elements from all sequences into resulting sequence
  • CancellationTokenSource to notify producers that early termination is requested

Instead let’s change behavior to allow producer proceed with getting next element from source sequence only once previously “produced” element is consumed (which basically means merged sequence enumerator reached it). This is sort of two way synchronization between producer and consumer where consumer cannot proceed unless there are ready elements and producer cannot proceed unless previous element is consumed. It is kind of similar to blocking collection. However in this case there is only one consumer and consumers aren’t blocked because the queue is full but rather until previously enqueued element is dequeued.Let’s code the thing!

class Consumer<T>
{
    private readonly IEnumerable<IEnumerable<T>> _sources;

    private readonly object _sync = new object();

    private readonly Queue<Producer<T>> _globalQueue = new Queue<Producer<T>>();
    private readonly Queue<Producer<T>> _localQueue = new Queue<Producer<T>>();

    // No need to mark these fields as volatile as 
    // only consumer thread updates and reads them
    private bool _done;
    private int _count;

    private readonly IList<Exception> _exceptions = new List<Exception>();

    public Consumer(IEnumerable<IEnumerable<T>> sources)
    {
        _sources = sources;
    }

    // Merges sources in parallel
    public IEnumerable<T> Merge()
    {
        Start();
        while (true)
        {
            Wait();
            while (_localQueue.Count > 0)
            {
                // Use local queue to yield values
                var producer = _localQueue.Dequeue();
                // Get the actual value out of ready producer
                // while simultaneously allowing it to proceed 
                // and observe early termination request or get 
                // next value
                var next = producer.GetNext(!_done);
                // Using Notificaiton<T> from Rx to simplify
                // sequence processing
                if (next.Kind != NotificationKind.OnNext)
                {
                    _count--;
                    if (next.Kind == NotificationKind.OnError)
                    {
                        // Observed exception leads to early 
                        // termination however merge will allow 
                        // already running (not waiting) 
                        // producers to finish (potentially 
                        // bringing more exceptions to be 
                        // observed) while requesting waiters 
                        // to terminate
                        _done = true;
                        // Store observed exception to be further 
                        // thrown as part of the aggregate 
                        // exception
                        _exceptions.Add(next.Exception);
                    }
                }
                else
                {
                    yield return next.Value;
                }
            }
            // Loop until all producers finished
            if (_count == 0)
            {
                // Either normally
                if (_exceptions.Count == 0)
                    yield break;
                // Or exceptions were observed
                throw new AggregateException(_exceptions);
            }
        }
    }

    // Notifies consumer of ready elements
    public void Enqueue(Producer<T> producer)
    {
        // Notify consumer of a ready producer
        lock (_sync)
        {
            // Consumer will either observe non-empty 
            // global queue (if it is not waiting already) 
            // or pulse (if it is waiting)
            _globalQueue.Enqueue(producer);
            Monitor.Pulse(_sync);
        }
    }

    // Waits for ready elements
    private void Wait()
    {
        lock (_sync)
        {
            // As the only consumer is draining the global queue
            // once an empty queue is observed and consumer is 
            // notified the queue will be non-empty
            if (_globalQueue.Count == 0)
                Monitor.Wait(_sync);
            // Copy whatever available to local queue to further 
            // drain without bothering taking a lock
            while (_globalQueue.Count > 0)
                _localQueue.Enqueue(_globalQueue.Dequeue());
        }
    }

    // Starts producers
    private void Start()
    {
        try
        {
            foreach (var source in _sources)
            {
                var producer = new Producer<T>(this, source.GetEnumerator());
                Task.Factory.StartNew(producer.Enumerate);
                _count++;
            }
        }
        catch (Exception ex)
        {
            // If none of producers are started successfully 
            // just rethrow the exception
            if (_count == 0)
                throw;
            // Some producers started notify them of early 
            // termination
            _done = true;
            // Store observed exception to be further thrown as 
            // part of the aggregate exception
            _exceptions.Add(ex);
        }
    }
}

class Producer<T>
{
    private readonly object _sync = new object();
    private readonly Consumer<T> _consumer;

    private IEnumerator<T> _enum;

    private volatile Notification<T> _next;
    private volatile bool _cont = true;
    private volatile bool _awake;

    public Producer(Consumer<T> consumer, IEnumerator<T> enumerator)
    {
        _consumer = consumer;
        _enum = enumerator;
    }

    // Drains source sequence
    public void Enumerate()
    {
        try
        {
            // Loop always observes non-null value as 
            // it is initially non-null and everytime consumer 
            // awakes producer it is set non-null value
            while (_cont && _enum.MoveNext())
            {
                // Set continuation flag to null and wait to be 
                // awoken by the consumer
                _awake = false;
                // Notify consumer of a ready element
                _next = new Notification<T>.OnNext(_enum.Current);
                _consumer.Enqueue(this);

                lock (_sync)
                {
                    // If consumer was pretty quick and producer 
                    // missed the pulse it will observe the awake 
                    // flag and thus won't wait
                    if (!_awake)
                        Monitor.Wait(_sync);
                }
            }
            _next = new Notification<T>.OnCompleted();
        }
        catch (Exception ex)
        {
            _next = new Notification<T>.OnError(ex);
        }
        finally
        {
            _enum.Dispose();
            _enum = null;
        }
        // Notify consumer that the producer completed 
        _consumer.Enqueue(this);
    }

    // Awakes/notifies producer that it can proceed
    public Notification<T> GetNext(bool cont)
    {
        // Store ready element in local variable 
        // because once producer is unleashed below 
        // it may override ready yet not returned element
        var next = _next;

        lock (_sync)
        {
            // Set awake flag in case producer miss the pulse 
            _awake = true;
            _cont = cont;
            Monitor.Pulse(_sync);
        }

        return next;
    }
}

With Producer<T> and Consumer<T> in place Merge combinator can be implemented as:

static class EnumerableEx
{
    public static IEnumerable<T> Merge<T>(IEnumerable<IEnumerable<T>> sources)
    {
        Contract.Requires(sources != null);
        return new Consumer<T>(sources).Merge();
    }
}

At any given moment at most n elements will be in a global queue where n is the number of source sequences.

Great exercise!

Sunday, February 13, 2011

Longest consecutive elements sequence

A tiny detail that can be uncovered by looking at a problem on a different angle usually is the key to the solution. So many times I looked at a great API design or problem solution saying “how I didn’t see it”. But sometimes I do see. It was a problem of searching for the longest consecutive elements sequence within unsorted array of integers. For example, in {5, 7, 3, 4, 9, 10, 1, 15, 1, 3, 12, 2, 11} sought sequence is {1, 2, 3, 4, 5}.

The first thing that comes mind is to sort given array in O(n log n ) and look for longest consecutive elements sequence. Eh, we can do better than that.

Using bit vector indexed by numbers from original array may not be justified due incomparable solution space and original array size (although time complexity is O(n)).

Let’s look at the problem more closely. The problem may be reduced to problem of effective range manipulation. Disjoint-set structure or interval trees offer O(n log n) building time complexity. But they do not take into consideration fact that we are dealing with integers. Knowing range boundaries we can definitely say what numbers are in there (for example, [1..3] contains 1, 2, 3). O(1) time complexity for range manipulation operations with O(1) space complexity for each range are the things we are looking for. We can do this using two hash tables:

  • ‘low’ with range start numbers as keys and range end numbers as values
  • ‘high’ with range end numbers as keys and range start numbers as values

For example, for a range [x..y] the tables will hold low[x] = y and high[y] = x. The algorithm looks the following:

  • Scan original array and for each element:
    • If it already belongs to any range skip it
    • Otherwise create unit size range [i..i]
    • If there is a range next to the right [i+1.. y] merge with it to produce [i..y] range
    • If there is a range next to the left [x..i-1] merge with it as well (either [i..i] or merged [i..y] will be merged with [x..i-1])
  • Scan through created ranges to find out the longest

The only question left is how to check if an element already belongs to some range as we are keeping only range boundaries in hash tables. Any number is either processed previously and thus in one of the ranges or a new range created out of it (and potentially merged with others). Thus if a number previously seen it is in some range and we do not need to process it. So before scanning original array we can simply remove any duplicates in O(n) time and space.

class Solver
{
    public static Tuple<int, int> FindMaxRange(IEnumerable<int> seq)
    {
        // Generate all ranges and select maximum
        return EnumerateRanges(seq)
            .Aggregate((x, y) => Length(x) > Length(y) ? x : y);
    }

    public static IEnumerable<Tuple<int, int>> EnumerateRanges(IEnumerable<int> seq)
    {
        var low = new Dictionary<int, int>();
        var high = new Dictionary<int, int>();

        // Remove duplicates
        foreach (var val in seq.Distinct())
        {
            // Create unit size range
            low[val] = high[val] = val;
            // Merge [i..i] with [i+1..y]
            var endsWith = MergeWithNext(val, low, high, 1);
            // Merge [i..endsWith] with [x..i-1]
            MergeWithNext(endsWith, high, low, -1);
        }

        return low.Select(p => Tuple.Create(p.Key, p.Value));
    }

    static int MergeWithNext(int currStart, IDictionary<int, int> low, IDictionary<int, int> high, int sign)
    {
        var currEnd = low[currStart];
        var nextStart = currEnd + sign;
        if (low.ContainsKey(nextStart))
        {
            low[currStart] = low[nextStart];
            high[low[currStart]] = currStart;
            low.Remove(nextStart);
            high.Remove(currEnd);
        }
        return low[currStart];
    }

    static int Length(Tuple<int, int> t)
    {
        return t.Item2 - t.Item1;
    }
}
The solution has O(n) time and space complexity assuming hash table implementation has O(1) time and space complexity for each operation/element.

Wednesday, January 5, 2011

Parallel string matching

String matching is about searching for occurrence (first or all occurrences – it makes difference from parallelization point of view as you shall see soon) of a pattern in a given text.

This problem naturally fits into data parallelism scenario although details depend on whether we want to find first or all occurrences of a pattern. Parallel searching for all occurrences is simpler as searching space is static (whole text needs to be looked through) in contrast to searching for the first occurrence that may be anywhere in a text (reducing search space improves performance otherwise the solution is either equal to sequential search till first occurrence or parallel search of all occurrences with getting the first one).

In order to find all pattern occurrences text can be separated into chunks that are processed in parallel. Each chunk overlaps with its immediate neighbor (except for the last one) for no more than length(pattern) - 1 characters to cover the case when pattern occurs in text such that chunk starting position is within pattern.

public static class ParallelAlgorithms
{
    // Find all occurrences of a pattern in a text
    public static IEnumerable<int> IndexesOf(this string text, string pattern, int startIndex, int count)
    {
        var proc = Environment.ProcessorCount;
        // Do range partitioning
        var chunk = (count + proc - 1) / proc;
        var indexes = new IEnumerable<int>[proc];

        Parallel.For(0, proc, 
            p => 
            {
                // Define overlapping with its immediate 
                // neighbor chunk except for the last one
                var before = p * chunk;
                var now = Math.Min(chunk + pattern.Length - 1, count - p * chunk);
                // Sequentially search for patterns occurences 
                // and store local result
                indexes[p] = SequentialIndexesOf(text, pattern, startIndex + before, now);
            });
        return indexes.SelectMany(p => p);
    }

    static IEnumerable<int> SequentialIndexesOf(string text, string pattern, int startIndex, int count)
    {
        for (var i = 0; i <= count - pattern.Length;)
        {
            var found = text.IndexOf(pattern, startIndex + i, count - i);
            // No pattern occurrences is found till the end
            if (found == -1)
                break;
            yield return found;
            // Proceed with next to found position
            i = found + 1;
        }
    }
}

Now for the first occurrence scenario search space must be reduced minimizing amount of unnecessary work (it will be non-zero in most cases due speculative processing).

One way to do this is to separate text into chunks and process them in parallel such that if pattern is found in chunk(i) processing of any chunk(j) where j > i is cancelled. Assuming n is the length of a text and p equals to logical processor count in the worst case (when first occurrence of a pattern is at the end of the first chunk) amount of the unnecessary work will be (p - 1)*n/p. On the other hand range partitioning has poor performance when workload is unbalanced.

What we can do instead is dynamic partitioning where whole text is separated into groups of chunks of the same length within group and chunk length in consequent group is two times large than previous group had. Group size should be set to number of logical processors. Thus, if c denotes chunk of unit size text will be separated like c, c, c, …, cc, cc, cc, …, cccc, cccc, cccc,…  Now this sequence of chunks can be processed in parallel (respecting order) breaking at first found occurrence. Thus in worst case at most p*c amount of unnecessary work will be done.

This partitioning strategy is used in PLINQ and called chunk partitioning. Although text can be treated as a sequence of characters and chunk partitioning can be used out of the box but that is not what we want as otherwise we will process individual characters rather than text chunks. Instead we’ll produce sequence of chunks manually and using single item partitioner and Parallel.ForEach process them in parallel respecting order.

public static class ParallelAlgorithms
{
    // Find first occurrence of a pattern in a text
    public static int IndexOf(this string text, string pattern, int startIndex, int count)
    {
        var minChunkSize = pattern.Length << 5;
        var maxChunkSize = minChunkSize << 3;
        // Create sequence of chunks
        var chunks = PartitionRangeToChunks(startIndex, count, minChunkSize, maxChunkSize, Environment.ProcessorCount);
        // Process chunks in parallel respecting order
        var chunkPartitioner = SingleItemPartitioner.Create(chunks);
        var concurrentBag = new ConcurrentBag<int>();
        Parallel.ForEach(chunkPartitioner,
            (range, loop) =>
            {
                var start = range.Item1;
                var length = Math.Min(startIndex + count - start, range.Item2 + pattern.Length - 1);
                var index = text.IndexOf(pattern, start, length);
                // No pattern occurrences in this chunk
                if (index < 0)
                    return;
                // Store shared result
                concurrentBag.Add(index);
                // Let running parallel iterations complete and 
                // prevent starting new ones
                loop.Break();
            });
        // Pick first occurrence or indicate no occurrence
        return concurrentBag.Count > 0 ? concurrentBag.Min() : -1;
    }

    static IEnumerable<Tuple<int, int>> PartitionRangeToChunks(int start, int count, int minChunkSize, int maxChunkSize, int doubleAfter)
    {
        var end = start + count;
        var chunkSize = Math.Min(minChunkSize, count);
        while (start < end)
        {
            for (var i = 0; i < doubleAfter && start < end; i++)
            {
                var next = Math.Min(end, start + chunkSize);
                yield return Tuple.Create(start, next - start);
                start = next;
            }

            chunkSize = Math.Min(maxChunkSize, chunkSize * 2);
        }
    }
}

Search in parallel! =)

Monday, December 27, 2010

Parallel graph search

Graph is a concept used to describe relationships between things where vertices (or nodes) represent things and edges represent relationships. Graph is made up of two finite sets (vertices and edges) and denoted by G = (V, E). Nodes usually are denoted by non-negative integers (for example, graph with 3 nodes with have nodes denoted by 0, 1, 2) and edges are pairs of nodes (s, t) where s is a source and t is a target node.

Searching for a simple path (path with no repeated vertices) that connects two nodes is one of the basic problems that are solved using graph search. We need to visit each node once (meaning visited nodes tracking is on) to find out is there any path between two given nodes. Time complexity for graph search is O(V).

Let’s parallelize the thing to speed it up.

Graph representation is irrelevant for our case so let’s assume implementation of the following public API is in place:

public class Graph
{
    // Initializes new graph with given number of nodes and 
    // set of edges
    public Graph(int nodes, IEnumerable<Tuple<int, int>> edges);
    // Gets number of nodes in the graph
    public int Nodes { get; }
    // Gets list of nodes adjacent to given one
    public IEnumerable<int> AdjacentTo(int node);
}

Starting with a source node for each node being processed search condition is checked and search is expanded to adjacent nodes that are not yet visited. Search continues until condition is met or no more reachable not visited nodes are left.

First part of the solution is to process nodes in parallel until no more is left or early termination is requested. The idea is to maintain two work queues for current and next phases. Items from current phase queue are processed in parallel and new work items are added to next phase queue. At the end of phase queues are switched. While this mechanism covers “while not empty” scenario, cooperative cancellation covers early search termination.

public static class ParallelAlgorithms
{
    public static void DoWhileNotEmpty<T>(IEnumerable<T> seed, Action<T, Action<T>> body, CancellationToken cancellationToken)
    {
        // Maintain two work queues for current and next phases
        var curr = new ConcurrentQueue<T>(seed);
        var next = new ConcurrentQueue<T>();
        var parallelOptions = new ParallelOptions {CancellationToken = cancellationToken};
        // Until there is something to do
        while (!curr.IsEmpty)
        {
            // Function to add work for the next phase
            Action<T> add = next.Enqueue;
            // Process current work queue in parallel while 
            // populating next one
            Parallel.ForEach(curr, parallelOptions, t => body(t, add));
            // Switch queues and empty next one
            curr = next;
            next = new ConcurrentQueue<T>();
        }
    }
}

Now that we have processing mechanism we need to bring visited nodes tracking mechanism. As part of our solution we must find a path that connects two given nodes if one exists. Thus instead of just remembering which nodes were visited we’ll keep track of parent nodes used to visit a node. Unvisited nodes are marked with –1. Once search is finished we can reconstruct path backwards between target and source by following parent links.

public static class ParallelGraph
{
    public static int[] Search(this Graph graph, int source, Func<Tuple<int, int>, bool> func)
    {
        const int undef = -1;
        // Start with dummy loop edge
        var seed = Enumerable.Repeat(Tuple.Create(source, source), 1);
        // Mark all nodes as not yet visited
        var tree = Enumerable.Repeat(undef, graph.Nodes).ToArray();
        tree[source] = source;
        // Cancellation token source used to exit graph search 
        // once search condition is met.
        var cancellationTokenSource = new CancellationTokenSource();
        try
        {
            // Until there are reachable not visited nodes 
            ParallelAlgorithms
                .DoWhileNotEmpty(seed, (edge, add) =>
            {
                var from = edge.Item2;

                // and search condition is not met
                if (!func(edge))
                    cancellationTokenSource.Cancel();

                // Try expand search to adjacent nodes
                foreach (var to in graph.AdjacentTo(from))
                {
                    // Expand search to not yet visited nodes 
                    // (otherwise if node is expanded and 
                    // then checked during processing we may 
                    // end up with lots of already visited 
                    // nodes in memory which is a problem 
                    // for large dense graphs) marking nodes 
                    // to visit next
                    if (Interlocked.CompareExchange(ref tree[to], from, undef) != undef)
                        continue;
                    add(Tuple.Create(from, to));
                }
            }, cancellationTokenSource.Token);
        }
        catch (OperationCanceledException ex)
        {
            // Check if exception originated from parallel 
            // processing cancellation request
            if (ex.CancellationToken != cancellationTokenSource.Token)
                throw;
        }
        // In the end we have spanning tree that connects
        // all the nodes reachable from source
        return tree;
    }

    public static IEnumerable<int> FindPath(this Graph graph, int from, int to)
    {
        // Search within graph until reached 'to'
        var tree = Search(graph, from, edge => edge.Item2 != to);
        // Check if there is a path that connects 'from' and 'to'
        if (tree[to] == -1)
            return Enumerable.Empty<int>();
        // Reconstruct path that connects 'from' and 'to'
        var path = new List<int>();
        while (tree[to] != to)
        {
            path.Add(to);
            to = tree[to];
        }
        path.Add(to);
        // Reverse path as it is reconstructed backwards
        return Enumerable.Reverse(path);
    }
}

Done!

Thursday, November 11, 2010

Joining Microsoft

Find a job you like and you add five days to every week.

- H. Jackson Brown, Jr.

I decided to join Microsoft in attempt to make my week longer. Lot’s of things to learn and challenging problems to solve. But this is where interesting topics come from.

I will continue chasing state of the art soon. Now I need to get materialized in Seattle area =).

Sunday, October 17, 2010

Parallel merge sort

Divide and conquer algorithm solves the problem by:

  1. dividing problem into two or more smaller independent sub-problems of the same type of problem
  2. solving each sub-problem recursively
  3. and combining their results.

Sub-problem independency makes divide and conquer algorithms natural for dynamic task parallelism. Quicksort is a good example (actually you can find its parallel version here). We’ll focus on parallel merge sort as it better parallelizes. Basically the term parallelism means ratio T1/Ti where T1 is the execution time on a single processor and Ti is the execution on infinite number of processors.

Let’s implement sequential merge sort first as parallel merge sort will use it. In general merge sort in terms of divide and conquer defined as follows:

  1. If the array is of size 0 or 1 it is already sorted; otherwise divide it into two sub-arrays of about half the size
  2. recursively sort each sub-array
  3. and merge sorted sub-arrays into one sorted array.

It is quite straightforward with a small improvement that allows to avoid unnecessary allocations and data copying.

class MergeSortHelper<T>
{
    private readonly IComparer<T> _comparer;

    public MergeSortHelper()
        : this(Comparer<T>.Default)
    {
    }

    public MergeSortHelper(IComparer<T> comparer)
    {
        _comparer = comparer;
    }

    public void MergeSort(T[] array, int low, int high, bool parallel)
    {
        // Create a copy of the original array. Switching between
        // original array and its copy will allow to avoid 
        // additional array allocations and data copying.
        var copy = (T[]) array.Clone();
        if (parallel)
            ParallelMergeSort(array, copy, low, high, GetMaxDepth());
        else
            SequentialMergeSort(array, copy, low, high);
    }

    private void SequentialMergeSort(T[] to, T[] temp, int low, int high)
    {
        if (low >= high)
            return;
        var mid = (low + high) / 2;
        // On the way down the recursion tree both arrays have 
        // the same data so we can switch them. Sort two 
        // sub-arrays first so that they are placed into the temp 
        // array.
        SequentialMergeSort(temp, to, low, mid);
        SequentialMergeSort(temp, to, mid + 1, high);
        // Once temp array contains two sorted sub-arrays
        // they are merged into target array.
        SequentialMerge(to, temp, low, mid, mid + 1, high, low);
        // On the way up either we are done as the target array
        // is the original array and now contains required 
        // sub-array sorted or it is the temp array from previous 
        // step and contains smaller sub-array that will be 
        // merged into the target array from previous step 
        // (which is the temp array of this step and so we 
        // can destroy its contents).
    }

    // Although sub-arrays being merged in sequential version 
    // are adjacent that is not the case for parallel version 
    // and thus sub-arrays boundaries must be specified 
    // explicitly.
    private void SequentialMerge(T[] to, T[] temp, int lowX, int highX, int lowY, int highY, int lowTo)
    {
        var highTo = lowTo + highX - lowX + highY - lowY + 1;
        for (; lowTo <= highTo; lowTo++)
        {
            if (lowX > highX)
                to[lowTo] = temp[lowY++];
            else if (lowY > highY)
                to[lowTo] = temp[lowX++];
            else
                to[lowTo] = Less(temp[lowX], temp[lowY])
                                ? temp[lowX++]
                                : temp[lowY++];
        }
    }

    private bool Less(T x, T y)
    {
        return _comparer.Compare(x, y) < 0;
    }
    
    ...

Now we need to parallelize it. Let’s get obvious out of the way. Sorting two sub-arrays can be done in parallel just like parallel quicksort. We can proceed to merge step only once both sub-arrays are sorted.

Parallelizing merge is a more interesting task. Sequential version looks like a single indivisible task. We need to find a way to separate merging into independent tasks that can be run in parallel. The idea behind algorithm is the following:

  1. Let’s assume we want to merge sorted arrays X and Y. Select X[m] median element in X. Elements in X[ .. m-1] are less than or equal to X[m]. Using binary search find index k of the first element in Y greater than X[m]. Thus Y[ .. k-1] are less than or equal to X[m] as well. Elements in X[m+1 .. ] are greater than or equal to X[m] and Y[k .. ] are greater. So merge(X, Y) can be defined as concat(merge(X[ .. m–1], Y[ .. k–1]), X[m], merge(X[m+1 .. ], Y[k .. ]))
  2. now we can recursively in parallel do merge(X[ .. m-1], Y[ .. k–1]) and merge(X[m+1 .. ], Y[k .. ])
  3. and then concat results.

Although algorithm description pretty abstract there will be more detailed comments in code below.

    ...

    private const int SEQUENTIAL_THRESHOLD = 2048;

    // Recursion depth is utilized to limit number of spawned 
    // tasks. 
    private void ParallelMergeSort(T[] to, T[] temp, int low, int high, int depth)
    {
        if (high - low + 1 <= SEQUENTIAL_THRESHOLD || depth <= 0)
        {
            // Resort to sequential algorithm if either 
            // recursion depth limit is reached or sub-problem 
            // size is not big enough to solve it in parallel.
            SequentialMergeSort(to, temp, low, high);
            return;
        }

        var mid = (low + high) / 2;
        // The same target/temp arrays switching technique 
        // as in sequential version applies in parallel 
        // version. sub-arrays are independent and thus can 
        // be sorted in parallel.
        depth--;
        Parallel.Invoke(
            () => ParallelMergeSort(temp, to, low, mid, depth),
            () => ParallelMergeSort(temp, to, mid + 1, high, depth)
            );
        // Once both taks ran to completion merge sorted 
        // sub-arrays in parallel.
        ParallelMerge(to, temp, low, mid, mid + 1, high, low, depth);
    }

    // As parallel merge is itself recursive the same mechanism
    // for tasks number limititation is used (recursion depth).
    private void ParallelMerge(T[] to, T[] temp, int lowX, int highX, int lowY, int highY, int lowTo, int depth)
    {
        var lengthX = highX - lowX + 1;
        var lengthY = highY - lowY + 1;

        if (lengthX + lengthY <= SEQUENTIAL_THRESHOLD || depth <= 0)
        {
            // Resort to sequential algorithm in case of small 
            // sub-problem or deep recursion.
            SequentialMerge(to, temp, lowX, highX, lowY, highY, lowTo);
            return;
        }

        if (lengthX < lengthY)
        {
            // Make sure that X range no less than Y range and 
            // if needed swap them.
            ParallelMerge(to, temp, lowY, highY, lowX, highX, lowTo, depth);
            return;
        }

        // Get median of the X sub-array. As X sub-array is 
        // sorted it means that X[lowX .. midX - 1] are less 
        // than or equal to median and X[midx + 1 .. highX] 
        // are greater or equal to median.
        var midX = (lowX + highX) / 2;
        // Find element in the Y sub-array that is strictly 
        // greater than X[midX]. Again as Y sub-array is 
        // sorted Y[lowY .. midY - 1] are less than or equal 
        // to X[midX] and Y[midY .. highY] are greater than 
        // X[midX].
        var midY = BinarySearch(temp, lowY, highY, temp[midX]);
        // Now we can compute final position in the target 
        // array of median of the X sub-array.
        var midTo = lowTo + midX - lowX + midY - lowY;
        to[midTo] = temp[midX];
        // The rest is to merge X[lowX .. midX - 1] with 
        // Y[lowY .. midY - 1] and X[midx + 1 .. highX] 
        // with Y[midY .. highY] preceeding and following 
        // median respectively in the target array. As 
        // pairs are idependent from their final position 
        // perspective they can be merged in parallel.
        depth--;
        Parallel.Invoke(
            () => ParallelMerge(to, temp, lowX, midX - 1, lowY, midY - 1, lowTo, depth),
            () => ParallelMerge(to, temp, midX + 1, highX, midY, highY, midTo + 1, depth)
            );
    }

    // Searches for index the first element in low to high 
    // range that is strictly greater than provided value 
    // and all elements within specified range are smaller 
    // or equal than index of the element next to range is 
    // returned.
    private int BinarySearch(T[] from, int low, int high, T lessThanOrEqualTo)
    {
        high = Math.Max(low, high + 1);
        while (low < high)
        {
            var mid = (low + high) / 2;
            if (Less(from[mid], lessThanOrEqualTo))
                low = mid + 1;
            else
                high = mid;
        }
        return low;
    }

    private static int GetMaxDepth()
    {
        // Although at each step we split unsorted array
        // into two equal size sub-arrays sorting them 
        // not be perfectly balanced because parallel merge 
        // may not be balanced. So we add some extra space for 
        // task creation and so will keep CPUs busy.
        return (int) Math.Log(Environment.ProcessorCount, 2) + 4;
    }
}

Now that the thing is ready let’s do some comparison. We’ll use small benchmark helper.

public class Benchmark
{
    public static IEnumerable<long> Run<T>(Func<int, T[]> generator, Action<T[]> action, int times)
    {
        var samples = new long[times];
        var sw = new Stopwatch();
        for (var i = 0; i < times; i++)
        {
            var input = generator(i);
            sw.Restart();
            action(input);
            sw.Stop();
            samples[i] = sw.ElapsedMilliseconds;
        }
        return samples;
    }
}

We'll calculate and compare average running time for a number of samples for sequential quicksort, parallel quicksort and parallel merge sort.

const int count = 10 * 1000 * 1000;
const int iterations = 10;
var rnd = new Random();

var input = new int[count];
Func<int, int[]> gen = i =>
                {
                    // Avoid allocating large array for every
                    // run and just fill existing one 
                    for (var j = 0; j < input.Length; j++)
                        input[j] = rnd.Next();

                    return input;
                };

Action<int[]> seqArraySort, parallelQuickSort, parallelMergeSort;
// If Array.Sort<T>(T[] a, IComparer<T> c) sees 
// Comparer<T>.Default it resorts to unmanaged CLR provided 
// sorting implementation, so we use here dummy comparer to 
// force it to use managed quicksort.
seqArraySort = a => Array.Sort(a, new IntComparer());
// Parallel quicksort implementation from TPL extras.
// We use dummy comparer as well because at some point
// parallel quicksort resorts back to Array.Sort
parallelQuickSort = a => ParallelAlgorithms.Sort(a, new IntComparer());
// Our parallel merge sort implementation
parallelMergeSort = a => new MergeSortHelper<int>().MergeSort(a, 0, a.Length - 1, true);

var sq = Benchmark.Run(gen, seqArraySort, iterations).Average();
var pq = Benchmark.Run(gen, parallelQuickSort, iterations).Average();
var pm = Benchmark.Run(gen, parallelMergeSort, iterations).Average();

On a two cores machine I got that parallel merge sort is more than 2x faster than sequential quicksort and up to 25% faster than parallel quicksort but at the cost of additional O(n) space. Still it is a good example of how to use dynamic task parallelism. Enjoy!

Sunday, September 19, 2010

Traverse binary tree in level order by spiral

Another puzzle is at stake, folks. This time it is binary tree related (not necessarily binary search tree as we are not interested in data relations but rather in binary tree structure). We need to traverse binary tree level by level (level is defined as set of all nodes at the same distance from root) in such a way that traversal direction within level changes from level to level thus forming a spiral. For example, consider binary tree below (it is rotated 90 degrees counter clockwise). Asterisk means NIL node.

                *
            20
                *
         19
                   *
               18
                   *
            17
                   *
               13
                   *
      12
          *
   11
       *
10
             *
          7
             *
       6
             *
          5
             *
    4
          *
       2
          *

Desired traversal for this tree will be: 10, 4, 11, 12, 6, 2, 5, 7, 19, 20, 17, 13, 18. Dissected into levels it looks like:

  1. 10
  2. 4, 11
  3. 12, 6, 2
  4. 5, 7, 19
  5. 20, 17
  6. 13, 18

It may not be obvious how to approach the problem. We’ll start with a well know problem of traversing binary tree level by level (also known as breadth first traversal). It can be done using queue.

static IEnumerable<T> LevelOrderLeftToRight<T>(TreeNode<T> root)
{
    if (root == null)
       yield break;

    var next = new Queue<TreeNode<T>>();
    next.Enqueue(root);

    while(next.Count > 0)
    {
       var node = next.Dequeue();
       yield return node.Data;

       EnqueueIfNotNull(next, node.Left);
       EnqueueIfNotNull(next, node.Right);
    }
}

static void EnqueueIfNotNull<T>(Queue<T> queue, T value)
   where T:class
{
    if (value != null)
       queue.Enqueue(value);
}

The queue contains yet to be examined nodes in level by level left to right (as we first enqueue left child and then right) order. At each step node at the front is dequeued, examined and its children are enqueued for further examination thus preserving level by level left to right order. It yields the following traversal: 10, 4, 11, 2, 6, 12, 5, 7, 19, 17, 20, 13, 18. Let’s split node sequence into levels:

  1. 10
  2. 4, 11
  3. 2, 6, 12
  4. 5, 7, 19
  5. 17, 20
  6. 13, 18

That looks close to sought-for order except every odd numbered level has must be reversed. With all nodes in the same container (queue) it is hard to do this. Let’s change code a little bit to make every two adjacent layers are separated into different containers.

static IEnumerable<T> LevelOrderLeftToRight<T>(TreeNode<T> root)
{
    if (root == null)
       yield break;

    var curr = new Queue<TreeNode<T>>();
    var next = new Queue<TreeNode<T>>();
    next.Enqueue(root);

   do
   {
      // Swap level containers
      Swap(ref curr, ref next);
      // Examine all nodes at current level
      while (curr.Count > 0)
      {
         var node = curr.Dequeue();
         yield return node.Data;
         // Fill next level preserving order
         EnqueueIfNotNull(next, node.Left);
         EnqueueIfNotNull(next, node.Right);
      }
      // Continue until next level has nodes
   } while (next.Count > 0);
}

static void Swap<T>(ref T a, ref T b)
{
   var tmp = a;
   a = b;
   b = tmp;
}

With adjacent levels separated we can change order within levels. Next level is a set of child nodes from current level. FIFO container (queue) makes sure that child nodes of earlier examined node will also be examined earlier. But this is opposite to what we are looking for. Change it to LIFO container (stack)! And that’s it. Child nodes of earlier examined node will be examined later.

static IEnumerable<T> LevelOrderBySpiral<T>(TreeNode<T> root)
{
   if (root == null)
      yield break;

   var curr = new Stack<TreeNode<T>>();
   var next = new Stack<TreeNode<T>>();
   // Specifies direction for the next level
   var leftToRight = true;
   next.Push(root);

   do
   {
      Swap(ref curr, ref next);
      while (curr.Count > 0)
      {
         var node = curr.Pop();
         yield return node.Data;
         // If next level must be traversed from left to right
         // we must first push right child node and then left
         // and in opposite order if next level will be 
         // traversed from right to left
         PushIfNotNull(next, leftToRight ? node.Right : node.Left);
         PushIfNotNull(next, leftToRight ? node.Left : node.Right);
      }
      // Change direction within level
      leftToRight = !leftToRight;
   } while (next.Count > 0);
}

static void PushIfNotNull<T>(Stack<T> stack, T value)
   where T : class
{
   if (value != null)
      stack.Push(value);
}

The code yields sought-for order.

Tuesday, August 24, 2010

External merge sort

If data to be sorted doesn’t fit into main memory external sorting is applicable. External merge sort can be separated into two phases:

  1. Create initial runs (run is a sequence of records that are in correct relative order or in other words sorted chunk of original file).
  2. Merge created runs into single sorted file.

To implement this algorithm I will use solutions from my previous posts so it may be helpful for you to look at them:

Let’s assume that M records at the same time are allowed to be loaded into main memory. One of the ways to create initial runs is to successively read M records from original file, sort them in memory and write back to disk. However we will use approach that allows us to create longer runs. It is called replacement selection.

The core structure behind this algorithm is priority queue. Taking one by one current minimum element out of the queue forms ordered sequence. And this is exactly what run stands for. The algorithm can be described as follows:

  1. Create two priority queues (that will contain items for current and next runs respectively) with capacity of M records.
  2. Prefill current run priority queue from unsorted file with M records.
  3. Create current run if there are elements in the current run priority queue:
    1. Take minimum element out of current run priority queue and append it to current run (basically write it to disk).
    2. Take next element from unsorted file (this is the replacement part of the algorithm) and compare it with just appended element.
    3. If it is less then it cannot be part of the current run (or otherwise order will be destroyed) and thus it is queued to the next  run priority queue.
    4. Otherwise it is part of the current run and it is queued to the current run priority queue.
    5. Continue steps 1 through 4 until current run priority queue is empty.
  4. Switch current and next runs priority queues and repeat step 3.

At any given moment at most M records are loaded into main memory as single written element into current run is replaced with single element from unsorted file if any (depending on comparison it either goes into current or next run).

Next step is to merge created initial runs. For the merge step we will use simplified algorithm (more advanced algorithms work with multiple physical devices to distribute runs, take into account data locality, etc.) based on k-way merge:

  1. Append created runs into a queue.
  2. Until there are more than one run in the queue:
    1. Dequeue and merge K runs into a single run and put it into the queue.
  3. Remaining run represents sorted original file.

Yeap, it is that simple. And let’s code it.

The implementation abstracts file structure and reading/writing details making algorithm more concise and easier to understand.

abstract class ExternalSorter<T>
{
	private readonly IComparer<T> m_comparer;
	private readonly int m_capacity;
	private readonly int m_mergeCount;

	protected ExternalSorter(IComparer<T> comparer, int capacity, int mergeCount)
	{
		m_comparer = comparer;
		m_capacity = capacity;
		m_mergeCount = mergeCount;
	}

	// Sorts unsorted file and returns sorted file name
	public string Sort(string unsorted)
	{
		var runs = Distribute(unsorted);
		return Merge(runs);
	}

	// Write run to disk and return created file name
	protected abstract string Write(IEnumerable<T> run);
	// Read run from file with given name
	protected abstract IEnumerable<T> Read(string name);

	// Merge step in this implementation is simpler than 
	// the one used in polyphase merge sort - it doesn't
	// take into account distribution over devices
	private string Merge(IEnumerable<string> runs)
	{
		var queue = new Queue<string>(runs);
		var runsToMerge = new List<string>(m_mergeCount);
		// Until single run is left do merge
		while (queue.Count > 1)
		{
			// Priority queue must not contain records more than 
			// required
			var count = m_mergeCount;
			while (queue.Count > 0 && count-- > 0)
				runsToMerge.Add(queue.Dequeue());
			// Perform n-way merge on selected runs where n is 
			// equal to number of physical devices with 
			// distributed runs but in our case we do not take 
			// into account them and thus n is equal to capacity
			var merged = runsToMerge.Select(Read).OrderedMerge(m_comparer);
			queue.Enqueue(Write(merged));

			runsToMerge.Clear();
		}
		// Last run represents source file sorted
		return queue.Dequeue();
	}

	// Distributes unsorted file into several sorted chunks
	// called runs (run is a sequence of records that are 
	// in correct relative order)
	private IEnumerable<string> Distribute(string unsorted)
	{
		var source = Read(unsorted);
		using (var enumerator = source.GetEnumerator())
		{
			var curr = new PriorityQueue<T>(m_comparer);
			var next = new PriorityQueue<T>(m_comparer);
			// Prefill priority queue to capacity which is used 
			// to create runs
			while (curr.Count < m_capacity && enumerator.MoveNext())
				curr.Enqueue(enumerator.Current);
			// Until unsorted source and priority queues are 
			// exhausted
			while (curr.Count > 0)
			{
				// Create next run and write it to disk
				var sorted = CreateRun(enumerator, curr, next);
				var run = Write(sorted);

				yield return run;

				Swap(ref curr, ref next);
			}
		}
	}

	private IEnumerable<T> CreateRun(IEnumerator<T> enumerator, PriorityQueue<T> curr, PriorityQueue<T> next)
	{
		while (curr.Count > 0)
		{
			var min = curr.Dequeue();
			yield return min;
			// Trying to move run to an end enumerator will 
			// result in returning false and thus current 
			// queue will simply be emptied step by step
			if (!enumerator.MoveNext())
				continue;

			// Check if current run can be extended with 
			// next element from unsorted source
			if (m_comparer.Compare(enumerator.Current, min) < 0)
			{
				// As current element is less than min in 
				// current run it may as well be less than 
				// elements that are already in the current 
				// run and thus from this element goes into 
				// next run
				next.Enqueue(enumerator.Current);
			}
			else
			{
				// Extend current run
				curr.Enqueue(enumerator.Current);
			}
		}
	}

	private static void Swap<U>(ref U a, ref U b)
	{
		var tmp = a;
		a = b;
		b = tmp;
	}
}

In the example below I created type that sorts text files containing single number per line.

class TextFileOfNumbersExternalSorter : ExternalSorter<int>
{
	public TextFileOfNumbersExternalSorter(int capacity, int mergeCount)
		: base(Comparer<int>.Default, capacity, mergeCount)
	{
	}

	protected override string Write(IEnumerable<int> run)
	{
		var file = Path.GetTempFileName();
		using (var writer = new StreamWriter(file))
		{
			run.Run(writer.WriteLine);
		}
		return file;
	}

	protected override IEnumerable<int> Read(string name)
	{
		using (var reader = new StreamReader(name))
		{
			while (!reader.EndOfStream)
				yield return Int32.Parse(reader.ReadLine());
		}
		File.Delete(name);
	}
}

That is used like this:

// capacity, mergeCount and unsortedFileName are initialized elsewhere
var sorter = new TextFileOfNumbersExternalSorter(capacity, mergeCount);
var sortedFileName = sorter.Sort(unsortedFileName);

That’s it folks!

Saturday, August 21, 2010

Minimum window that contains all characters

Like most of life's problems, this one can be solved with bending

- Bender B.Rodrigues

Let’s bend another problem. Given set of characters P and string T find minimum window in T that contains all characters in P. Applicable solution is restricted to O(length(T)) time complexity. For example, given a string T “of characters and as” and set of characters T in a form of a string “aa s” the minimum window will be “and as”.

The problem can be broken into two parts:

  • How to select window?
  • How to check that selected window contains all characters from P?

Selecting every possible window (all unique pairs (i, j) where 0 <= i <= j < length(T)) will lead to solution worse than O(length(T)^2) because you still need to check if all characters from P are within selected window. Instead we will check every possible window ending position. Thus there are at least length(T) windows to consider.

Any feasible window has length equal to or greater than length(P). Performing recheck for any considered window will result in a solution no better than O(length(T)*length(P)). Instead we need to use check results from previous iteration.

Now we need to make sure that checking if a particular character is in P is done in an optimal way. Taking into account that a particular character may appear more than once and window thus must contain appropriate number of characters. We will use hash table to map unique characters from P to their count for fast lookup.

And now let’s tie all things together.

  • Until reached the end of T move by one current window ending position.
  • Append next character to the end of previous window which to this moment doesn’t contain all necessary characters. Char to count map is used to track the number of characters left to find. Basically if character is in P count is decremented. The number may become negative meaning that there are more than required characters.
  • If unmatched character count goes to zero the window contains all required characters. However there may be redundant characters. Thus we try to compact current window. It is ok to do this as we are looking for minimum window and any window that is extended from this one won’t be better.
  • Once window is compacted compare it with the minimum one and updated it if needed.
  • If current window contains all the characters remove from it the first one simply by moving by one starting position to make sure that at each iteration previous window doesn’t contain all the characters (there is no point in appending new characters to a window that already contains all of the required ones).

Code the thing! =)

static string FindMinWindow(string t, string p)
{
	// Create char to count mapping for fast lookup
	// as some characters may appear more than once
	var charToCount = new Dictionary<char, int>();
	foreach (var c in p)
	{
		if (!charToCount.ContainsKey(c))
			charToCount.Add(c, 0);
		charToCount[c]++;
	}

	var unmatchesCount = p.Length;
	int minWindowLength = t.Length + 1, minWindowStart = -1;
	int currWindowStart = 0, currWindowEnd = 0;
	for (; currWindowEnd < t.Length; currWindowEnd++)
	{
		var c = t[currWindowEnd];
		// Skip chars that are not in P
		if (!charToCount.ContainsKey(c))
			continue;
		// Reduce unmatched characters count
		charToCount[c]--;
		if (charToCount[c] >= 0)
			// But do this only while count is positive
			// as count may go negative which means 
			// that there are more than required characters
			unmatchesCount--;

		// No complete match, so continue searching
		if (unmatchesCount > 0)
			continue;

		// Decrease window as much as possible by removing 
		// either chars that are not in T or those that 
		// are in T but there are too many of them
		c = t[currWindowStart];
		var contains = charToCount.ContainsKey(c);
		while (!contains || charToCount[c] < 0)
		{
			if (contains)
				// Return character to P
				charToCount[c]++;

			c = t[++currWindowStart];
			contains = charToCount.ContainsKey(c);
		}

		if (minWindowLength > currWindowEnd - currWindowStart + 1)
		{
			minWindowLength = currWindowEnd - currWindowStart + 1;
			minWindowStart = currWindowStart;
		}

		// Remove last char from window - it is definitely in a 
		// window because we stopped at it during decrease phase
		charToCount[c]++;
		unmatchesCount++;
		currWindowStart++;
	}

	return minWindowStart > -1 ?
	       t.Substring(minWindowStart, minWindowLength) :
	       String.Empty;
}

Every character is examined at most twice (during appending to the end and during compaction) so the whole solution has O(length(T)) time complexity assuming hash table lookup is O(1) operation.

Tuesday, August 17, 2010

Print numbers by spiral

Recently I came across simple yet interesting coding problem. So here is the deal. You are given positive integer N. Print first N ^ 2 positive integers in matrix form in a such a way that within matrix numbers form spiral starting from its center and goring clockwise. For example, for N = 5 matrix to be printed is:

21 22 23 24 25
20 7 8 9 10
19 6 1 2 11
18 5 4 3 12
17 16 15 14 13

Optimize it for speed and space.

One way you can approach it is to create N x N matrix and fill it with numbers that form spiral and then print whole matrix row by row. But this solution will be of N ^ 2 space complexity. Let’s try to reach O(1) space complexity.

The key observation here is how matrix changes when N changes by 1.

N = 1.

1

N = 2.

1 2
4 3

N = 3.

7 8 9
6 1 2
5 4 3

N = 4.

7 8 9 10
6 1 2 11
5 4 3 12
16 15 14 13

Can you see the pattern here? At every step we extend previous matrix (P) with additional column and row (C). If N is even we extend previous matrix of size N – 1 with right column and bottom row

P C
C C

and with left column and top row if it is odd

C C
C P

This leads us to naturally recursive algorithm. We have three cases:

  1. Print whole row of the current matrix (top when N is odd or bottom when N is even).
  2. Print row from previous matrix of size N - 1 first and then print value that belongs to current matrix (when N is even).
  3. Print value that belongs to current matrix and then print row from previous matrix of size N - 1 (when N is odd).
  4. Print matrix line by line.

So basically to print a row we need to know matrix size N and row index. Here goes the solution.

static void Print(int n)
{
	for(int i = 0; i < n; i++)
	{
		PrintLine(n, i);
		Console.WriteLine();
	}
}

static void PrintLine(int n, int i)
{
	// Number of integers in current matrix
	var n2 = n*n;
	// Number of itegers in previous matrix of size n - 1
	var m2 = (n - 1)*(n - 1);

	if (n % 2 == 0)
	{
		if (i == n - 1)
		{
			// n is even and we are at the last row so just 
			// print it
			for(int k = n2; k > n2 - n; k--)
			{
				PrintNum(k);
			}
		}
		else
		{
			// Print row from previous matrix of size n - 1 
			// first and then print value that belongs to current 
			// matrix. Previous matrix is at the top left corner 
			// so no need to adjust row index
			PrintLine(n - 1, i);
			// Skip all integers from previous matrix and upper 
			// ones in this columnas integers must form clockwise 
			// spiral
			PrintNum(m2 + 1 + i);
		}
	}
	else
	{
		if (i == 0)
		{
			// n is odd and we are at the first row so just 
			// print it
			for(int k = m2 + n; k <= n2; k++)
			{
				PrintNum(k);
			}
		}
		else
		{
			// Print value that belongs to current matrix and
			// then print row from previous matrix of size n - 1
			// Skip all integers from previous matric and bottom
			// ones in this column as integers must form clockwise
			// spiral
			PrintNum(m2 + n - i);
			// Previous matrix is at the bottom right corner so
			// row index must be reduced by 1
			PrintLine(n - 1, i - 1);
		}
	}
}

static void PrintNum(int n)
{
	Console.Write("{0, -4}  ", n);
}

If stack is not considered then this solution has O(1) space complexity otherwise O(N).