Showing posts with label programming. Show all posts
Showing posts with label programming. Show all posts

Tuesday, January 31, 2012

Hadoop 1.0

There is finally a 1.0 version of hadoop. One of my biggest complaints using hadoop since version 0.10 is that something breaks with every release and the usual retort is that compatibility will come with 1.0, and it seemed like 1.0 was just around the corner for years. I haven't been using hadoop as much for the last six months, so I didn't notice that there was finally a 1.0 release until today even though it was announced towards the end of last year. As a user this news inspires some hope that there might be a hadoop upgrade that just works. We'll have to see as new versions come out.

That said, it is a little dissapointing to see some of the cruft that is still around in the APIs. In particular org.apache.hadoop.mapred
and org.apache.hadoop.mapreduce packages are both still around. So if I want to write a map reduce job which API should I use? Even worse the javadocs still don't provide much clarity on which of these APIs should be preferred and neither appears to be deprecated. Maybe they have good reasons, but this adds a lot of confusion for users and, in my selfish opinion, should have been cleaned up before a 1.0 release.

The other big problem I've had with hadoop is trying to keep a stack of hadoop, pig, oozie, hbase, etc all working together flawlessly and being able to update individual components without too much worry on whether the rest of the stack will play nice. This is much easier to do if hadoop provides clean, simple, and well documented APIs that these other tools can build on. At first glance, the 1.0 apidocs look like they just slapped a 1.0 label on the current pile of crud and did not remove any of the old garbage and deprecated APIs.

If they actually maintain compatibility between 1.x releases it will still be a win for users, but hopefully for 2.0 they focus on a clean simple API for users and get rid of a bunch of old cruft. It would also be nice if we don't have to wait 6 years for 2.0.

Tuesday, November 29, 2011

Damn auto-completion

I really wish people would be more careful when using the auto-completion features of IDEs. Though auto-completion does provide some time savings it is also a frequent contributor to sloppy coding. After updating some packages our system broke because someone had incorrectly imported com.google.inject.internal.Sets instead of com.google.common.collect.Sets. Is it too much to ask for people to at least look at the code that is getting generated?

Saturday, November 12, 2011

Top-K Selection

Ok, so I need to select the top-K values from a list of size N where K is much smaller than N. Two approaches that immediately come to mind are:

  1. Sort the list and select the first K elements. Running time O(N*lg(N)).
  2. Insert the elements of the list into a heap and pop the top K elements. Running time O(N + K*lg(N)).

As a variant on option 2 a colleague proposed using a small heap that would have at most K elements at any given time. When the heap is full an element would be popped off before a new element could be added. So if I wanted the top-K smallest values I would use a max-heap and if the next value from the input is smaller than the largest value on the heap I would pop off the largest value and insert the smaller value. The running time for this approach is O(N*lg(K)).

For my use case both N and K are fairly small. The size of N will be approximately 10,000 elements and K will typically be 10, but can be set anywhere from 1 to 100. A C++ test program that implements all three approaches and tests them for various sizes of K is provided at the bottom of this post. The table below shows the results for K equal to 10, 50, and 100. You can see that all three approaches have roughly the same running time and increasing the size of K has little impact on the actual running time.

K=10K=50K=100
sort0.4689770.4669180.467177
fullheap0.0993250.1028390.106735
smallheap0.0180630.0289480.040435

Here is the chart for all values of K:


So clearly the third approach with the small heap is the winner. With a small K it is essentially linear time and an order of magnitude faster than the naive sort. The graph shows both the average time and the 99th-percentile so you can see the variation in times is fairly small. This first test covers the sizes for my use case, but out of curiosity, I also tested the more interesting case with a fixed size K and varying the size of N. The graph for N from 100,000 to 1,000,000 in increments of 100,000 tells the whole story:


Source code:

#include <algorithm>
#include <cstdlib>
#include <ctime>
#include <iostream>
#include <numeric>
#include <vector>

using namespace std;

typedef void (*topk_func)(vector<int> &, vector<int> &, int);

bool greater_than(const int &v1, const int &v2) {
    return v1 > v2;
}

void topk_sort(vector<int> &data, vector<int> &top, int k) {
    sort(data.begin(), data.end());
    vector<int>::iterator i = data.begin();
    int j = 0;
    for (; j < k; ++i, ++j) {
        top.push_back(*i);
    }
}

void topk_fullheap(vector<int> &data, vector<int> &top, int k) {
    make_heap(data.begin(), data.end(), greater_than);
    for (int j = 0; j < k; ++j) {
        top.push_back(*data.begin());
        pop_heap(data.begin(), data.end(), greater_than);
        data.pop_back();
    }
}

void topk_smallheap(vector<int> &data, vector<int> &top, int k) {
    for (vector<int>::iterator i = data.begin(); i != data.end(); ++i) {
        if (top.size() < k) {
            top.push_back(*i);
            push_heap(top.begin(), top.end());
        } else if (*i < *top.begin()) {
            pop_heap(top.begin(), top.end());
            top.pop_back();
            top.push_back(*i);
            push_heap(top.begin(), top.end());
        }
    }
    sort_heap(top.begin(), top.end());
}

void run_test(const string &name, topk_func f, int trials, int n, int k) {
    vector<double> times;

    for (int t = 0; t < trials; ++t) {
        vector<int> data;
        for (int i = 0; i < n; ++i) {
            data.push_back(rand());
        }

        clock_t start = clock();
        vector<int> top;
        f(data, top, k);
        clock_t end = clock();
        times.push_back(1000.0 * (end - start) / CLOCKS_PER_SEC);
    }

    sort(times.begin(), times.end());
    double sum = accumulate(times.begin(), times.end(), 0.0);
    double avg = sum / times.size();
    double pct90 = times[static_cast<int>(trials * 0.90)];
    double pct99 = times[static_cast<int>(trials * 0.99)];
    cout << name << " "
         << k << " "
         << avg << " "
         << pct90 << " "
         << pct99 << endl;
}

int main(int argc, char **argv) {

    cout << "method k avg 90-%tile 99-%tile" << endl;

    int trials = 1000;
    int n = 10000;
    for (int k = 1; k <= 100; ++k) {
        run_test("sort", topk_sort, trials, n, k);
        run_test("fullheap", topk_fullheap, trials, n, k);
        run_test("smallheap", topk_smallheap, trials, n, k);
    }

    return 0;
}

Friday, November 11, 2011

Broken pipe

If you are using python to write a script, please properly handle the broken pipe exception. We have set of command line tools at work that are extremely annoying to use because if you pipe the output through other standard tools, e.g., head, it spits out a worthless exception about a broken pipe. Consider the following test script:

#!/usr/bin/env python
import sys
buffer = ""
for i in range(100000):
    buffer += "%d\n" % i
sys.stdout.write(buffer)

Pipe it into head and look at the output:

$ ./test.py | head -n1
0
ERROR:root:damn
Traceback (most recent call last):
  File "./test.py", line 6, in <module>
    sys.stdout.write(buffer)
IOError: [Errno 32] Broken pipe

I know there was a broken pipe and I don't care. Just swallow the worthless exception so I can see the meaningful output. This is probably the number one reason I often mutter "damn python crap" when using some of these tools. So if you are writing scripts in python, please be considerate and handle the broken pipe exception. Here is an example for quick reference:

#!/usr/bin/env python
import errno
import sys
try:
    buffer = ""
    for i in range(100000):
        buffer += "%d\n" % i
    sys.stdout.write(buffer)
except IOError, e:
    if e.errno != errno.EPIPE:
        raise e

Scala, regex, and null

One of the perks of using scala is that I hardly ever see a NullPointerException unless I'm working with java libraries. The primary reason is because most scala libraries tend to use Option rather than null. However, while using a regex with pattern matching I was surprised by a NullPointerException when trying to look at the result of an optional capture group. Consider the following example:

scala> val Pattern = "(a)(b)?".r
Pattern: scala.util.matching.Regex = (a)(b)?

scala> "a" match { case Pattern(a, b) => printf("[%s][%s]%n", a, b) }
[a][null]

scala> "ab" match { case Pattern(a, b) => printf("[%s][%s]%n", a, b) }
[a][b]

I just assumed that b would be of type Option[String]. There is probably a good reason for this travesty, my guess would be something about making it work with the type system, but after using scala for a while it just seems wrong to be getting a null value.

Monday, September 26, 2011

Man in the middle

You typically hear the expression man in the middle in the context of an attack where someone is actively eavesdropping on communications that are intended to be private. However, it is also an invaluable tool for debugging network programs. One of my favorite tools is netcat and this tool makes it trivial to implement a simple eavesdropping script. This can also be done with tools such as tcpdump, but I find that netcat is a bit simpler for most tasks and it is more likely to be available on the machine.

The script I typically use is shown below. Essentially it has an endless loop that listens on a port, tees the data that comes in to a request log file, sends the input to a remote server, tees the response from the remote server to a response log, and then writes the data back to a named pipe that is connected to stdin on the netcat process that was listening.

#!/bin/bash

function logFile {
   echo "$(date +%Y-%m-%d-%H-%M-%S).${1}.log"
}

function serveRequests {
   port=$1
   remoteHost=$2
   remotePort=$3
   while true; do
       rm -f backpipe
       mkfifo backpipe
       cat backpipe |
           nc -l $port |
           tee -a $(logFile request) |
           nc $remoteHost $remotePort |
           tee -a $(logFile response) >backpipe
   done
}

port=${1:-12345}
remoteHost=${2:-localhost}
remotePort=${3:-80}
serveRequests $port $remoteHost $remotePort

Wednesday, September 21, 2011

Inexcusable laziness

Below is the text of a warning email I received from an internal intrusion detection system:
Subject: IMPORTANT: security violations found for cluster

Security violations found for instances of cluster: foobar

To see the full report go to:
http://ids.mycompany.com/reports/cluster/<clustername>
There are many things I think could be improved with this email, but the primary thing that annoyed me was the link to the full report. Why not insert the actual cluster name so I can just click on the link? The way it is I have to copy the prefix of the url and then type or copy in the cluster it is complaining about. Are functional links too much to ask for?

Monday, February 21, 2011

zfcat

While messing with java, I have frequently wanted an easy way to cat a file inside of a zip archive such as a jar, war, ear, and whatever other names they have cooked up for a zip file with a manifest. As one example, when trying to setup maven to generate an OSGi bundle, I want to generate the jar file and then look at the generated manifest to see if it is correct. There is probably an existing command line tool that provides this functionality, but I couldn't find one with the little bit of searching I tried and it is trivial to write such a tool. So I wrote a quick "zip file cat" tool called zfcat in Scala:
#!/bin/sh
exec scala $0 $@
!#

import java.io._
import java.util.zip._

object zfcat {
    def copy(in: InputStream, out: OutputStream): Unit = {
        val buffer = new Array[Byte](4096)
        var length = in.read(buffer)
        while (length != -1) {
            out.write(buffer, 0, length)
            length = in.read(buffer)
        }
    }

    def main(args: Array[String]): Unit = {
        if (args.length < 2) {
            System.err.printf("Usage: zfcat <archive> [file ...]%n")
            exit(1)
        }

        val zfile = new ZipFile(new File(args(0)))
        args.tail.foreach(file => {
            val entry = zfile.getEntry(file)
            if (entry != null) {
                val in = zfile.getInputStream(entry)
                copy(in, System.out)
                in.close
            } else {
                System.err.printf("Warning: zfcat: %s: no such file%n", file)
            }
        })
    }
}

zfcat.main(args)
This tool worked ok, but it just felt really slow compared to other command line tools. Java in general, and languages that run off of the JVM such as Scala, seems to be a terrible choice for writing command line tools. The problem is that the JVM is designed for long running tasks. For command line tools where most of the time it will be a very short lived job, programs based on the JVM are just too slow. So I wrote a new version of zfcat in Go:
package main

import (
    "archive/zip"
    "fmt"
    "io"
    "os"
)

func die(err os.Error) {
    if err != nil {
        fmt.Fprintf(os.Stderr, "Error: %s\n", err)
        os.Exit(1)
    }
}

func copy(r io.Reader, w io.Writer) os.Error {
    const BUFFER_SIZE = 4096
    var buffer [BUFFER_SIZE]byte
    for {
        switch nr, er := r.Read(buffer[:]); true {
            case nr <  0: return er
            case nr == 0: return nil
            case nr >  0: if nw, ew := w.Write(buffer[0:nr]); nw != nr {
                return ew
            }
        }
    }
    return nil
}

func main() {
    if len(os.Args) < 3 {
        fmt.Fprintf(os.Stderr, "Usage: %s <archive> [file ...]\n", os.Args[0])
        os.Exit(1)
    }

    reader, err := zip.OpenReader(os.Args[1])
    die(err)

    var zfiles = make(map[string] *zip.File)
    for i := range reader.File {
        file := reader.File[i]
        zfiles[file.FileHeader.Name] = file
    }

    files := os.Args[2:]
    for i := range files {
        file := zfiles[files[i]]
        if file != nil {
            r, err := file.Open()
            die(err)
            die(copy(r, os.Stdout))
            die(r.Close())
        } else {
            fmt.Fprintf(os.Stderr, "Warning: %s: %s: no such file\n",
                os.Args[0], files[i])
        }
    }
}
The Go version is fast for short lived jobs and seems to work great. I compared the times for three versions 1) running Scala as a script, 2) Scala pre-compiled, and 3) compiled Go version. The times to cat a small manifest file:
$ time ./zfcat.scala test.jar META-INF/MANIFEST.MF
Manifest-Version: 1.0
Created-By: 1.6.0_22 (Apple Inc.)


real    0m1.588s
user    0m1.076s
sys     0m0.083s
$ time scala zfcat test.jar META-INF/MANIFEST.MF
Manifest-Version: 1.0
Created-By: 1.6.0_22 (Apple Inc.)


real    0m0.635s
user    0m0.822s
sys     0m0.070s
$ time ./zfcat test.jar META-INF/MANIFEST.MF
Manifest-Version: 1.0
Created-By: 1.6.0_22 (Apple Inc.)


real    0m0.028s
user    0m0.003s
sys     0m0.018s
Pre-compiling more than doubles the speed of the Scala version, but it still takes over half a second. The Go version takes less than 30ms. I'm not sure I like the use of the error return values in Go and having to explicitly check them all the time. The decision to limit the use of exceptions as a control structure was deliberate, but I haven't bothered to read through the arguments for their proposal and it is a minor annoyance right now. I was also a little disappointed that I couldn't run valgrind, apparently the 6g compiler generates an executable that valgrind doesn't understand. Valgrind might work with executables generated by gccgo, but I didn't try it.
$ valgrind ./zfcat test.jar META-INF/MANIFEST.MF
bad executable (__PAGEZERO is not 4 GB)
valgrind: ./zfcat: cannot execute binary file

Saturday, February 5, 2011

Quicksort

I recently came across C.A.R. Hoare's 1962 paper on Quicksort. The paper is divided into two parts: theory and implementation. Most of the theory portion will be familiar to anyone that has taken an algorithms course. There are a few differences in terminology, the most noticeable to me was that Hoare uses the term bound for what is usually called the pivot in modern presentations. I also learned a new phrase, mutatis mutandis, that according to wikipedia means:
Mutatis mutandis is a Latin phrase meaning "by changing those things which need to be changed" or more simply "the necessary changes having been made".
What is really interesting is the implementation details and what they reveal about the state of computers at the time. The first example that comes up is the use of an alternative sort when a partition gets small enough. This is still a common quicksort optimization, but what is interesting is the size at which he recommends this optimization:
Furthermore, if a segment consists of less than three or four items (depending on the characteristics of the computer), then it will be advantageous to sort it by the use of a program specially written for sorting a particular small number of items.
I remember having to compare insertion sort with quick sort in college to find an appropriate cutoff point and as I recall it was around 50 items. Of course, this will vary with the computer, nuances of the sort implementations, and the data set being tested. I just would have expected this number to be at least 10 or so. Another example, is the data provided in table 1:

Number of ItemsMerge SortQuicksort
5002 min 8 sec1 min 21 sec
1,0004 min 48 sec3 min 8 sec
1,5008 min 15 sec*5 min 6 sec
2,00011 min 0 sec*6 min 47 sec

The asterisk by the results for the last two rows of merge sort indicate that those figures had to be estimated because they could not be computed. The reason provided is the limited store size available on the computer:
These figures were computed by formula, since they cannot be achieved on the 405 owing to limited store size.
The description of the machine and its capabilities:
The National-Elliott 405 computer has a delay-line working store of 512 locations, and a magnetic-disc backing store of 16,384 words. The average access time for the working store is 0.8 msec and the average access time for a block of 64 words in the backing store is 32 msec. There are 19 words of immediate-access storage, which are used to contain instructions and working space of the inner loops: the time taken by such loops is about 0.15 msec per instruction.
To get a better idea of the kind of computer he is talking about, see this video from the Yorkshire film archive showing the installation of a 405 computer at Dansom Lane. It won't help for understanding the memory limitations, but it is still interesting to get an idea of what the physical machine they were testing on looked like. The paper is not too clear on exactly what the data set looked like other than the single sentence:
The figures relate to six-word items with a single-word key.
Based on this description it is clear that the problem with merge sort is that it does not sort in-place. Duplicating the data set would mean that merge sort on the 405 could support a maximum data set size of 16,384 words / 2 copies / 6 words per item = 1,365 items. Unfortunately, there is no description of the formula used to calculate the last two times for merge sort. I was curious what the times would be for something like insertion sort. Insertion sort is in-place so it should be able to handle the larger data sets, but what would the expected running time be? An extremely rough approximation would be to assume exactly N2 instructions. If we use the number from the description in the paper of 0.15 msec per instruction and assume that for each instruction we would need to have at least one read from the working store at a cost of 0.8 msec, then we can use a cost of 0.95 msec per operation. This is almost certainly an underestimate of the actual time it would take. If nothing else I'm completely ignoring access to the magnetic disc. Crunch the numbers and look what we wind up with:

Number of ItemsN lg(N)N2
5000 min 4 sec3 min 57 sec
1,0000 min 9 sec15 min 50 sec
1,5000 min 15 sec35 min 37 sec
2,0000 min 20 sec1 hour 3 min 20 sec

Over an hour to sort 2000 items compared to less than 7 minutes for quicksort.

Friday, January 14, 2011

Damn TCLLIB htmlparse

I recently needed to write a quick script to extract some information from HTML documents. Poking around I found the machine already had TCLLIB installed so I thought it would be a good opportunity to try out the htmlparse library. The library is extremely easy to use except for one annoying pitfall, the handling of attributes. The parser will invoke a callback function with the tag name, text, and attributes. A default function is provided called ::htmlparse::debugCallback. From the documentation of the param argument:
The fourth argument, param, contains the un-interpreted list of parameters to the tag.
I thought for sure they must be kidding. Do they really expect me to parse the attribute data myself? The reason I'm using a library is so I don't have to worry about all the intricacies of HTML parsing. I decided to give it a try with the example below:
#!/usr/bin/env tclsh

package require htmlparse 1.2

set html {
<html>
  <head><title>Test HTML Page</title></head>
  <body>
    <p>This is some test <a target = "_blank"
    href="http://w3.org/html"         >HTML</a> content.</p>
  </body>
</html>
}

::htmlparse::parse $html
The output of running this example:
$ ./htmlparse.tcl 
==> hmstart {} {} {
}
==> html {} {} {
  }
==> head {} {} {}
==> title {} {} {Test HTML Page}
==> title / {} {}
==> head / {} {
  }
==> body {} {} {
    }
==> p {} {} {This is some test }
==> a {} {target = "_blank"
    href="http://w3.org/html"         } HTML
==> a / {} { content.}
==> p / {} {
  }
==> body / {} {
}
==> html / {} {
}
==> hmstart / {} {}
Sure enough, the attributes are all in one big string just as the documentation stated. This is one of those times I was hoping the documentation was wrong.

Wednesday, January 5, 2011

Integer representation of dates

For a system at work, I needed to represent dates and times as integers. An integer representation is required because we need to import the data into a search index so we can efficiently perform point and range queries. The search system only supports range queries for integer types and one of our main use cases is querying based on a date range. Further, the system does not have a specific time data type, it just supports signed integers of 1, 2, 4, and 8 bytes. So the problem is how should we represent a date or time as an integer?

First, lets make the requirements more explicit. The input data comes is a JSON feed with ISO 8601 formatted strings for date and timestamp fields. Whatever representation we choose should be able to represent the full range of dates from the input including %Y, %Y-%m, %Y-%m-%d, and %Y-%m-%dT%H:%M:%S (formats specified using the notation for strftime). In addition, it is preferable if the representation is easy to use for a person, i.e., I would like to be able to encode a known date without needing a program.

Problems with Unix Time

The obvious first question is can we just use Unix time? It is well known and widely supported, but there are a few problems trying to use it for our application.

The first problem is that it has a limited range. If we use a 32-bit integer we run into the year 2038 problem. Right now we don't really care about dates that far in the future, but we do care about historical dates so the year 1901 problem is an immediate concern. This problem is easily solved by using a 64-bit integer giving a range of 1970 plus or minus around 290 billion years. We can afford to be that short-sighted.

The second problem is for some fields we do not know the precise timestamp. We only know a particular year, month, or day. Somehow we need to represent these dates as being distinct from a precise date. For example if all we know is that something occurred in 1970 we don't want to encode it as 1970-01-01, because now we can't distinguish it from something that we knew happened on January 1, 1970. Unix time uses all of the bits already so there is no room to represent these vague dates.

A third problem is that it is quite cumbersome for a person to figure out the unix timestamps for a given range. Say I want to find all entries within the year 1980, I would need to compute the unix timestamp values for the beginning (1980-01-01T00:00:00 is 315532800) and end of the year (1980-12-31T23:59:59 is 347155199). This is trivial to code, but for a person doing it manually it is a pain.

The second problem is the killer, Unix time is not an option. The other issues are annoying, but we could live with them, especially considering the benefits of using a representation that is common and widely supported like unix time. So, what is the alternative?

Concatenated Fields Representation

This is a representation I've seen numerous times, but I wasn't able to find an official name or spec for it. The idea is pretty simple, we have a date represented as %Y-%m-%dT%H:%M:%S, we just remove the punctuation symbols and use the values as digits so we get the numeric value %Y%m%d%H%M%S. The first thing that becomes apparent is that we need an integer type that can support at least 14 decimal digits (4 year + 2 month of year + 2 day of month + 2 hour of day + 2 minute of hour + 2 second of minute = 14 digits). If we only care about the date we need at least 8 digits and for just the time we need 6 digits. A 16-bit signed integer as a range of -32768 to 32767 so we need at least a 32-bit integer. The table below shows the max value for 32 and 64 bit signed integers using this representation for a time, date, and datetime.

32-bit max64-bit max
Time23:59:5923:59:59
Date214748-12-31922337203685477-12-31
DateTimetoo small922337203-12-31T23:59:59

For a time or date a 32-bit signed integer is more than enough. Time in particular is limited so the max value will be 23:59:59 and using 64-bits provides no value. For date and datetime types the additional room can be used to expand the range of years supported. However, even with 32-bit integers the date type can works through the year 214748 which is more than enough for our use cases, but if needed 64-bits will get you past the year 922 trillion. For datetime a 32-bit integer is too small for 14 decimal digits. A 64-bit value will provide support past the year 922 million, and more importantly it covers the full range of the 32-bit date.

How does this fit our requirements? It covers the full range of time values that we expect to get in the input. Vague dates that are missing the day and/or month can be encoded as %Y0000. The zeros are not valid for the month or the day of month so we can clearly distinguish from more specific dates. This format is also pretty easy for a human to use. So it fixes the problems with unix time, but have we added any new problems?

Yes. Probably the most annoying problem is that not all values are valid. For example the date 19701332 is nonsense. There is no 13th month or 32nd day of a month. This representation makes it really easy to have values that are complete garbage. In fact, most of the distinct values for the integer representation will not be valid dates. The table below shows the percentage of the distinct values that are valid:

32-bit max64-bit max
Time0.002% (86,400 of 4,294,967,296)0.000000000000468% (86,400 of 18,446,744,073,709,551,616)
Date2.08% (89,335,584 of 4,294,967,296)2.08% (383,692,276,733,158,848 of 18,446,744,073,709,551,616)
DateTimetoo small12.00% (2,213,609,288,860,213,248 of 18,446,744,073,709,551,616)

The second problem is closely related to the first. This representation wastes a lot of space and doesn't make good use of all the bits. For my purposes this is not a big problem, but if you are looking for a compact representation then this is a bad choice.

Finally, we cannot make use of the negative values. In the ISO 8601 representation the year can be negative to specify BCE dates. This approach cannot be used with this representation because the sign would apply to the entire value, not just the year, and that breaks the ordering. For example -0004-01-04 is earlier than -0004-12-30, but when represented as an integer -40104 is greater than -41230.

Summary

The concatenated fields representation is a user friendly way to represent common era dates as an integer. It is particularly useful if you need to support dates where only partial information is present, e.g., just the year.

Sunday, November 28, 2010

Automatic Generation of Color Palettes

A problem I've had several times is how to automatically generate colors for use in graphs. The user will be able to select some set of items that should be included and I need to select colors for each item. The requirements are:
  • The color for an item should be easy to distinguish from other items. Of course, it would also be nice if the colors looked decent, but from a functional perspective, the requirement is to be able to distinguish the items in the graph.
  • I need to be able to generate an arbitrary number of colors and avoid a fixed color palette with a predefined number. For practical purposes the number will be limited by the ability to distinguish different colors, but it would be nice for the mechanism to scale gracefully as the number of items increases.
  • The color of the background, for my purposes white, cannot be used.
When I examined a few tools, I found that most worked by having a fixed set of colors. My first attempt was to perform a naive increment of RGB pixel values. A trivial increment works well for shades of gray. This produces a palette like:

Hex24816
000000                                                                
0E0E0E                                                                   
1C1C1C                                                                    
2A2A2A                                                                    
383838                                                                    
464646                                                                    
545454                                                                    
626262                                                                    
707070                                                                    
7E7E7E                                                                    
8C8C8C                                                                    
9A9A9A                                                                    
A8A8A8                                                                    
B6B6B6                                                                    
C4C4C4                                                                    
D2D2D2                                                                    

The problem is that grayscale can be difficult to distinguish with more than a few colors. That is why tools like gnuplot use line patterns and shapes. However, most of my use cases are for graphs shown on a color monitor so there is no need to limit to grayscale. What happens if we try a naive increment with color? My first attempt was to treat the color as a three byte integer and simply divide the desired number of colors to get the increment value. Looking at the palette below you can see the results are poor:

Hex24816
000000                                                                    
0FFFF0                                                                    
1FFFE0                                                                    
2FFFD0                                                                    
3FFFC0                                                                    
4FFFB0                                                                    
5FFFA0                                                                    
6FFF90                                                                    
7FFF80                                                                    
8FFF70                                                                    
9FFF60                                                                    
AFFF50                                                                    
BFFF40                                                                    
CFFF30                                                                    
DFFF20                                                                    
EFFF10                                                                    

After looking around for a bit I found that the HSV representation is fairly well suited for this problem. HSV stands for hue, saturation, and value. The color space is represented as a cylinder:

For more background the paper Color Spaces for Computer Graphics gives a good overview and discusses how the various color spaces were designed with respect to human perception of color. To generate a palette the saturation and value settings can be fixed. The 360o for the hue can be divided by the desired number of colors and then we just increment the angle for each color. This technique gives a nice palette, but for more than around 8 colors it will be difficult for a person to distinguish some shades.

Hex24816
FF0000                                                                    
FF5F00                                                                    
FFBF00                                                                    
DFFF00                                                                    
7FFF00                                                                    
1FFF00                                                                    
00FF3F                                                                    
00FF9F                                                                    
00FFFF                                                                    
009FFF                                                                    
003FFF                                                                    
1F00FF                                                                    
7F00FF                                                                    
DF00FF                                                                    
FF00BF                                                                    
FF005F                                                                    

The Scala code I used for generating the palettes is shown below.
object Colors {

    import java.awt.Color

    def grayscale(num: Int): Seq[Color] = {
        // Truncate the full range of values to make sure we can distinguish
        // from the color white, i.e., (256, 256, 256).
        val range = 256 - 32

        // Determine how much to increment for each color.
        val delta = range / num
        if (delta == 0) {
            throw new IllegalArgumentException(
                "grayscale can support at most " + range + " colors")
        }

        // Generate the sequence of colors
        (0 until num).map(n => {
            val value = n * delta
            new Color(value, value, value)
        })
    }

    def naiveIncrement(num: Int): Seq[Color] = {
        // Truncate the full range of values to make sure we can distinguish
        // from the color white, i.e., (256, 256, 256).
        val range = 0xFFFFFF - 0xFF

        // Determine how much to increment for each color.
        val delta = range / num
        if (delta == 0) {
            throw new IllegalArgumentException(
                "naive increment can support at most " + range + " colors")
        }

        // Generate the sequence of colors
        (0 until num).map(n => {
            val value = n * delta
            new Color((value >> 16) & 0xFF, (value >> 8) & 0xFF, value & 0xFF)
        })
    }

    def hsv(num: Int): Seq[Color] = {
        // Range is 360 degrees for the hue
        val range = 360.0

        // Determine how much to increment for each color.
        val delta = range / num
        if (delta < 1.0) {
            throw new IllegalArgumentException(
                "hsv can support at most " + range + " colors")
        }
        // Generate the sequence of colors
        (0 until num).map(n => {
            val hue = n * delta
            val h = hue / 60.0
            val x = ((1 - Math.abs(h % 2 - 1)) * 255).toInt
            val c = h match {
                case h if 0.0 <= h && h < 1.0 => (255, x, 0)
                case h if 1.0 <= h && h < 2.0 => (x, 255, 0)
                case h if 2.0 <= h && h < 3.0 => (0, 255, x)
                case h if 3.0 <= h && h < 4.0 => (0, x, 255)
                case h if 4.0 <= h && h < 5.0 => (x, 0, 255)
                case h if 5.0 <= h && h < 6.0 => (255, 0, x)
                case _ => (0, 0, 0)
            }
            new Color(c._1, c._2, c._3)
        })
    }

    def main(args: Array[String]): Unit = {
        if (args.length < 2) {
            println("Usage: scala Colors <palette> <num>")
            exit(1)
        }

        // Supported palettes
        val palettes = Map(
            "grayscale" -> grayscale _,
            "naive"     -> naiveIncrement _,
            "hsv"       -> hsv _
        )

        // Generate colors and print
        palettes(args(0))(args(1).toInt).foreach(c => {
            println(c.getRGB.toHexString.toUpperCase.substring(2))
        })
    }
}