Greedy algorithm for activity selection

A typical example of problem that has an optimal solution by implementing a greedy algorithm is the activity selection one, here is its description on wikipedia. In few words, we have a bunch of activities, identified by a start and end time, and we want to find a maximum selection of non-conflicting elements.

A couple of test cases (written in C++11 for GoogleTest) should clarify the problem:

typedef std::pair<int, int> Activity;
typedef std::vector<Activity> Activities;

TEST(ActSel, Simple)
{
  Activities input { {1, 2}, {5, 9}, {0, 6}, {8, 9}, {3, 4}, {5, 7} };

  Activities output = selectMax(input);
  ASSERT_EQ(4, output.size());
  for(unsigned i = 1; i < output.size(); ++i)
    ASSERT_LE(output[i-1].second, output[i].first);
}

TEST(ActSel, Simple2)
{
  Activities input { {1, 4}, {3, 5}, {0, 6}, {3, 9}, {5, 9}, {5, 7}, {6, 10}, {8, 11}, {8, 12}, {2, 14}, {12, 16} };

  Activities output = selectMax(input);
  ASSERT_EQ(4, output.size());
  for(unsigned i = 1; i < output.size(); ++i)
    ASSERT_LE(output[i-1].second, output[i].first);
}

In both cases I expect a selection of four Activity objects in output. In the first case these elements: (1, 2) (3, 4) (5, 7) (8, 9), in the second one (1, 4) (5, 7) (8, 12) (12, 16), or maybe (8, 11) instead of (8, 12). As you can see, there could be more solutions, and the problem doesn't require you to be particolary choosy. Once you maximize the number of selected items, the actual value of each of them is not an issue.

Still, I want to ensure in my test cases that I peak a valid solution, so I check, through ASSERT_LE, that all the elements in the extracted sequence are ordered as expected.

As said above, this problem has a greedy optimal solution. What we have to do is just sorting the input elements by their second component (the end time), and then greedily accepting all the elements we can. As in this implementation:

Activities selectMax(Activities& input) // 1
{
  std::sort(input.begin(), input.end(), [](Activity a, Activity b) { return a.second < b.second; }); // 2

  Activities output;
  output.push_back(input[0]); // 3

  for(unsigned i = 0, j = 1; j < input.size(); ++j) // 4
  {
    if(input[j].first >= input[i].second) // 5
    {
      output.push_back(input[j]);
      i = j;
    }
  }

  return output;
}

1. We don't mind if this function modify the input parameter, so it is passed as non-constant reference. Beware that this should be know and accepted by the callers.
2. The "normal" STL sort() function would order the passed sequence by its first component. So we need to use the overloaded version that let as pass a predicate to be used as comparator. Using a C++11 lambda function, as shown here, makes it simple and elegant.
3. The first element is always selected.
4. Now we are ready to loop on all the other elements in the sequence. The real looping variable is j, while i is used to keep track of the last accepted element.
5. The first element after the last accepted one that starts not before the end of it, is pushed in the output sequence.

Go to the full post

Rod cutting by dynamic programming

A typical problem that suits well to show how dynamic programming works. We have a rod sized up to, let's say, 10. We can freely cut it in pieces (integer sized) to sell them at the best price. Given a price table, find out the way to get the most from it.

Here is a C++11 test case for GoogleTest that should clarify the requirements:

typedef std::vector<int> Vector;

unsigned cutRod(const Vector& price, unsigned size);

TEST(CutRod, Simple)
{
  Vector price { 1, 5, 8, 9, 10, 17, 17, 20, 24, 30 };

  ASSERT_EQ(30, cutRod(price, 10));
  ASSERT_EQ(25, cutRod(price, 9));
  ASSERT_EQ(18, cutRod(price, 7));
  ASSERT_EQ(10, cutRod(price, 4));
}

Given that price list, we see immediately that if we have in input a rod sized up to 3, the best strategy is selling it in a single piece.
But if we have a rod sized four, selling it untouched we'll get 9. Better if we split it in two rodes both sized two, that give us 5 + 5 = 10.
Similarly, a rod sized 7 is priced 17. If we split it in two parts sized 6 and 1, we'll get 17 + 1 = 18.

Brute force

We may think to apply a recursive approach to this problem to check all the possible cut combinations we can think of. It is very easy to write the relative code, but can't we expect it to scale well:

unsigned cutRod(const Vector& price, unsigned size)
{
  unsigned result = 0;
  for(unsigned i = 0; i < size; ++i)
    result = std::max(result, price[i] + cutRod(price, size - (i+1)));

  return result;
}

It is just a matter of recursively calling our function reducing each time the size of the rod we are considering. We compare any time the partial result with the one we have previously stored, keeping just the best one.

Top-down dynamic programming

One obvious problem in the previous solution is that we solve again and again the same sub-problems. We could save lot of running time storing them in a buffer. This simple but effective idea is the basic of the dynamic programming technique.

In this context, the bargain of using space to avoid spending time repeating the same task to get a partial result is called memoization (as keeping a memo).

Here is a possible top-down implementation, very close to the naive version seen above:

unsigned cutRod(const Vector& price, unsigned size)
{
  Vector memo(size + 1, -1);
  memo[0] = 0;

  return memoCutRod(price, size, memo);
}

Here cutRod() just creates a memo vector that would store the values for each sub-problem, as soon as we get its result. Then it would start the recursion calling a support function.

Notice that the memo buffer has one element more than the price list. This is for storing also the value of the dummy cut sized zero. It is not a strict necessity, since we know that it won't cost anything, but it would help to make our code cleaner.

unsigned memoCutRod(const Vector& price, unsigned size, Vector& memo)
{
  if(memo[size] >= 0) // 1
    return memo[size];

  unsigned result = 0; // 2
  for(unsigned i = 0; i < size; ++i)
    result = std::max(result, price[i] + memoCutRod(price, size - (i+1), memo));

  return memo[size] = result; // 3
}

1. If the realtive memo buffer is not negative, we have already calculated it. Job already done.
2. Otherwise we calculate the best price as seen before.
3. And we set a memo before returning it.

Bottom-up approach

Again dynamic programming, still using memoization as we have just seen, but starting from the bottom of the problem and crawling up to its top. In this case the implementation is even simpler, and avoid us the pain and the cost of recursion:

unsigned cutRod(const Vector& price, unsigned size)
{
  Vector memo(size + 1); // 1
  for(unsigned i = 1; i <= size; ++i) // 2
  {
    int value = -1; // 3
    for(unsigned j = 0; j < i; ++j) // 4
      value = std::max(value, price[j] + memo[i-j-1]);
    memo[i] = value;
  }

  return memo.back(); // 5
}

1. As in the top-down approach, we get an extra element in the memo vector, just to keep simpler the code. But this time we don't need to initialize it to a "bad" values, because we are setting it up iteratively starting from the beginning.
2. First element in memo is already set to its expected value (that is, zero) as courtesy of the vector constructor. We need to calculate all the other elements, up to the rightmost one.
3. Initialize the current memo value to less than the minimum acceptable value (meaning, less than zero).
4. Basically it is the same loop we have seen in the previous implementations, but here we explicitly go for the smaller element first.
5. End of the story, the answer is stored in the rightmost memo element.

Check on github for full C++11 code.

Go to the full post

Quicksort

Quicksort is known to be a fast O(N lg N) divide and conquer sorting algorithm, in its average behavior. Still we have to pay attention to the worst case scenario, that brings it to a O(N ** 2) time cost.

The idea is repetitively partitioning the data collection, splitting it in two parts, in a way that a randomly chosen pivot would be equal or greater than the values on its left partition, and then call again the quicksorting procedure, until there is nothing more left to sort. As one could easily spot, is a possible bad choice of the pivot that could lead to poor performances.

The resulting code should be something like this:

void quicksort(std::vector<int>& data, int left, int right) // 1
{
  if(left < right) // 2
  {
    int pivot = partition(data, left, right); // 3
    quicksort(data, left, pivot - 1); // 4
    quicksort(data, pivot + 1, right);
  }
}

1. The function requires in input the collection on which it should operate and the indexes of its leftmost and rightmost elements.
2. Check if the interval is not empty.
3. Split the original interval in two parts. On the left side we have all the values less or equal to the value in the pivot element.
4. Call again quicksort on the left and right partitions. Notice that the pivot element is already in the right place, and don't need to be considered anymore.

We just need to partition a collection as expected:

int partition(std::vector<int>& data, int left, int right)
{
  int pivot = data[right]; // 1
  int index = left - 1; // 2

  for(int i = left; i < right; ++i) // 3
  {
    if(data[i] <= pivot) // 4
      std::swap(data[++index], data[i]);
  }

  std::swap(data[++index], data[right]); // 5
  return index;
}

1. OK, this doesn't look smart. As pivot we always select the rightmost element in the interval.
2. Initialize index to the first-before-beginning position in the interval.
3. Loop on all the items in the interval, but the last one (that is, the pivot).
4. If the current element value is less than the pivot, let's swap it with the first not already used element on the left.
5. Finally, we swap the pivot (rightmost value in the interval) with the element next to index.

Full C++ code on github.

Go to the full post

Heapsort

Heapsort is an in-place sorting algorithm, like insertion sort, that asympthotically scores a nice O(N lg N) time complexity, like merge sort.

It makes use of the heap data structure, that is a normal array seen as a nearly complete binary tree, in its max-heap flavor. Meaning that its biggest value is placed in the first element of the array (considered as the root of the tree).

Implementing heapsort in C++ is pretty trivial, since it just a matter of calling two STL algorithm functions:

#include <vector>
#include <algorithm>

void heapsort(std::vector<int>& data)
{
  std::make_heap(data.begin(), data.end()); // 1
  std::sort_heap(data.begin(), data.end()); // 2
}

1. This make_heap() call rearranges the passed elements as a max-heap.
2. This sort_heap() call assumes that the passed sequence is a max-heap and sort it in ascending order.

But let's have some fun reimplementing by hand these two functions:

typedef std::vector<int> Vector;

void heapsort(Vector& data)
{
  buildMaxHeap(data);
  sortHeap(data);
}

We'll need a way to navigate down the binary heap:

unsigned childLeft(unsigned i) { return (2 * i) + 1; }
unsigned childRight(unsigned i) { return (2 * i) + 2; }

The root element is at index 0. Its children are on 1 and 2.
The left child of the root (index 1) has its own children on 3 and 4; its sibling (index 2) on 5 and 6.
We can get the index of the children of a generic node in a binary heap just multiplying its index by two and adding 1 (for the left one) or 2 (for the right one).
And we'll need a function to ensure that a node in the data structure is complying to the binary max-heap requisite (it should be bigger than its children):

void maxHeapify(Vector& data, unsigned i, unsigned len) // 1
{
  unsigned left = childLeft(i);
  unsigned right = childRight(i);

  unsigned largest = (left < len && (data[left] > data[i])) ? left : i;
  if(right < len && (data[right] > data[largest]))
    largest = right;

  if(largest != i) // 2
  {
    std::swap(data[i], data[largest]); // 3
    maxHeapify(data, largest, len); // 4
  }
}

1. We pass to the function the data collection, the index of the element that we are checking, and the number of element in the heap.
2. We have compared the current node value against the ones of its left and right children. If the largest one is a children, the rule of the heap is currently violated. We need to rearrange the nodes.
3. Firstly, we need to swap the nodes so that the largest one is above the other ones.
4. Then, we need to ensure that the swapping has not corrupted the max-heap structure.

We are finally ready to implement the two main functions:

void buildMaxHeap(Vector& data) // 1
{
  for(int i = data.size() / 2; i >= 0; --i) // 2
    maxHeapify(data, i, data.size());
}

void sortHeap(Vector& heap) // 3
{
  for(int i = heap.size() - 1; i > 0; --i) // 4
  {
    std::swap(heap[0], heap[i]); // 5
    maxHeapify(heap, 0, i); // 6
  }
}

1. Given an arbitrary collection of values, convert it to a max-heap.
2. Start from the bottom, up to the root.
3. We assume that the passed data respect the max-heap constrains.
4. We scan the heap starting from the rightmost element up to the second one.
5. We know that the heap root is the biggest element in the collection, so we swap it to the rightmost position.
6. Before starting a new iteration, we ensure that the data collection (except the one we have already sorted) is still a max-heap.

Full C++ source code and a couple of test cases for Google test on github.

Go to the full post