Casting to Integers Considered Harmful August 6th, 2009
Patrick Stein

Background

Many years back, I wrote some ambient music generation code. The basic structure of the code is this: Take one queen and twenty or so drones in a thirty-two dimensional space. Give them each random positions and velocities. Limit the velocity and acceleration of the queen more than you limit the same for the drones. Now, select some point at random for the queen to target. Have the queen accelerate toward that target. Have the drones accelerate toward the queen. Use the average distance from the drones to the queens in the i-th dimension as the volume of the i-th note where the notes are logarithmically spaced across one octave. Clip negative volumes to zero. Every so often, or when the queen gets close to the target, give the queen a new target.

It makes for some interesting ambient noise that sounds a bit like movie space noises where the lumbering enemy battleship is looming in orbit as its center portion spins to create artificial gravity within.

I started working on an iPhone application based on this code. The original code was in C++. The conversion to Objective C was fairly straightforward and fairly painless (as I used the opportunity to try to correct my own faults by breaking things out into separate functions more often).

Visualization troubles

uniform
The original code though chose random positions and velocities from uniform distributions. The iPhone app is going to involve visualization as well as auralization. The picture at the right here is a plot of five thousand points with each coordinate selected from a uniform distribution with range [-20,+20]. Because each axis value is chosen independently, it looks very unnatural.

gauss
What to do? The obvious answer is to use Gaussian random variables instead of uniform ones. The picture at the right here is five thousand points with each coordinate selected from a Gaussian distribution with a standard-deviation of 10. As you can see, this is much more natural looking.

How did I generate the Gaussians?

I have usually used the Box-Muller method of generating two Gaussian-distributed random variables given two uniformly-distributed random variables:

(defun random-gaussian ()
  (let ((u1 (random 1.0))
        (u2 (random 1.0)))
    (let ((mag (sqrt (* -2.0 (log u1))))
          (ang (* 2.0 pi u2)))
      (values (* mag (cos ang))
              (* mag (sin ang))))))

But, I found an article online that shows a more numerically stable version:

(defun random-gaussian ()
  (flet ((pick-in-circle ()
           (loop as u1 = (random 1.0)
                as u2 = (random 1.0)
                as mag-squared = (+ (* u1 u1) (* u2 u2))
                when (< mag-squared 1.0)
                return (values u1 u2 mag-squared))))
    (multiple-value-bind (u1 u2 mag-squared) (pick-in-circle)
      (let ((ww (sqrt (/ (* -2.0 (log mag-squared)) mag-squared))))
        (values (* u1 ww)
                (* u2 ww))))))

For a quick sanity check, I thought, let’s just make sure it looks like a Gaussian. Here, I showed the code in Lisp, but the original code was in Objective-C. I figured, If I just change the function declaration, I can plop this into a short C program, run a few thousand trials into some histogram buckets, and see what I get.

The trouble with zero

So, here comes the problem with zero. I had the following main loop:

#define BUCKET_COUNT 33
#define STDDEV       8.0
#define ITERATIONS   100000

  for ( ii=0; ii < ITERATIONS; ++ii ) {
    int bb = val_to_bucket( STDDEV * gaussian() );
    if ( 0 <= bb && bb < BUCKET_COUNT ) {
      ++buckets[ bb ];
    }
  }

I now present you with three different implementations of the val_to_bucket() function.

int val_to_bucket( double _val ) {
  return (int)_val + ( BUCKET_COUNT / 2 );
}

int val_to_bucket( double _val ) {
  return (int)( _val + (int)( BUCKET_COUNT / 2 ) );
}

int val_to_bucket( double _val ) {
  return (int)( _val + (int)( BUCKET_COUNT / 2 ) + 1 ) - 1;
}

As you can probably guess, after years or reading trick questions, only the last one actually works as far as my main loop is concerned. Why? Every number between -1 and +1 becomes zero when you cast the double to an integer. That’s twice as big a range as any other integer gets. So, for the first implementation, the middle bucket has about twice as many things in it as it should. For the second implementation, the first bucket has more things in it than it should. For the final implementation, the non-existent bucket before the first one is the overloaded bucket. In the end, I used this implementation instead so that I wouldn’t even bias non-existent buckets:

int val_to_bucket( double _val ) {
  return (int)lround(_val) + ( BUCKET_COUNT / 2 );
}

90% Accuracy Fails July 27th, 2009
Patrick Stein

Bruce Schneier pointed out this BBC News article about the base rate fallacy. In that article, they consider a Terrorist Detector that is right 90% of the time. Given that the test signals a positive for someone, how sure are you the person in question is a terrorist?

The gist of the math

Well, the math in that article is confusing. They say the answer is about 0.3%. But, they later do a diagram that makes it look like 1 in 100,000 to me. Further, they tacitly assume that while 10% of the time the test will identify a non-terrorist as a terrorist, it will never miss a real terrorist.

Regardless, the goal of the article is to show Bayes’s Theorem in action for something that has a low prior probability. In layman’s terms, if something is really, really unlikely, then you are very bad off if your test is only really good. It’s really unlikely that any given person is a terrorist. As such, if your test thinks every tenth person is a terrorist, your test is almost never right when it thinks it has found a terrorist.

In the article, they ask for help describing this fallacy in ways that will drive the point home to people on a gut level. I can kind of see how one could miss the idea in cancer screening. But, with terrorism?

The real world

Remember the last time you were on a plane? [Okay, I do… it was last Thursday.] You were sitting there with about a hundred fellow passengers. How many of them do you think were terrorists? My guess is that you were pretty sure (even before the flight took off) that zero of them were terrorists. How many of them bombed your flight? I’m pretty sure that number was zero (and not just because you wouldn’t be reading this otherwise). How many of you did the test think were terrorists? About ten.

What proportion of hundred-passenger flights would the test think are terrorist-free? Fewer than one in every 35,000 flights. Buy a lottery ticket if you’re on a flight where this (hypothetical) test thinks there are no terrorists because it just may be your lucky day.

To do the full analysis, you have to consider how many of those 35,000 flights really do have terrorists on them. But, it should be blatantly obvious to everyone that the answer is way, way under 34,999.

What is a terrorist, anyway?

This brings up another question that I always have when people start talking about detecting terrorists or the number of terrorists in the U.S., etc. When is someone a terrorist? If Khalid Sheikh Mohammed were traveling to Hawaii for some R&R, is he a terrorist? would the test detect him? should he be allowed on the plane?

Detecting terrorists is often talked about like detecting diseases and such. But, really, it is about detecting intentions. Intentions are far more evanescent than microbes.

Numbers 30:2, for programmers July 20th, 2009
Patrick Stein

The role-playing game Paranoia was set in a dystopian future where a paranoid computer had seized control to protect us from Commie, mutant traitors. In the game, secret societies were strictly forbidden. As such, every character was a member of some secret society or other. One of those secret societies was the FCCCP: First Church of Christ, Computer Programmer. This, of course, tickled me. Heck, it still does.

Outside the world of role-playing games, there is a long tradition of analyzing sacred texts, pulling out every possible thread of meaning and applying those meanings to every aspect of life. I’ve been pondering a series of articles that combines the exegesis with the computer programming.

Well, this week’s sermon at our synagogue was about the verse Numbers 30:2:

If a man makes a vow to the LORD, or swears an oath to bind himself by some agreement, he shall not break his word; he shall do according to all that proceeds out of his mouth.

How could I sit there without pondering Programming by Contract? In fact, I would say that this verse goes much beyond just the contract. Every line of code makes promises. Are they promises the code can keep?

To ground the discussion, I dug around for a short function that I wrote more than a year ago. The following snippet of code takes two arrays: counter and maximum. The array elements of each are assumed to be integers. The arrays are assumed to be the same length. The value of the i-th entry in counter is assumed to be less than the i-th entry in maximum. If you think of an axis-aligned, n-dimensional cube with one corner at the origin and its opposite corner at the maximum, the goal is for counter to scan through every set of integer coordinates in the cube (well, including the faces that touch the origin, but not including the faces that touch the maximum). If every entry in maximum is b, then the result should be just like incrementing a base b number that is displayed least-significant digit first. If maximum were \left( 60, 60, 24, 7, 4 \right), the counter would count out the seconds in February (non-leap year).

(defun increment-counter (counter maximum)
  (unless (null counter)
    (let ((nn (length counter)))
      (let ((result (make-array nn :initial-contents counter)))
        (dotimes (ii nn nil)
          (when (< (incf (aref result ii)) (aref maximum ii))
            (return result))
          (setf (aref result ii) 0))))))

What promises does the above code make? The function name increment-counter is a pretty obvious promise. The function is going to increment some sort of counter. Given just the function name, what would we expect? I would expect that it takes a single argument, that argument is a reference to an integer, and that integer is incremented in place. I would be way off. Our promise is already on sketchy grounds. I’m not coming up with a particularly good name at the moment though: step-cube-walker? increment-n-dimensional-counter?

The function doesn’t increment the counter in place. Instead, it returns a newly allocated counter (or nil, if we’ve reached the end). The function declaration doesn’t make any promise that maximum will be left untouched, though Lisp doesn’t make it easy to declare that. Lisp usually handles that sort of thing by convention. The function declaration also makes no claim that the counter and maximum arguments should be arrays, let alone that they should be the same length. (Technically, the code only requires that maximum be an array. Counter can be any sequence. Further, maximum is allowed to be longer than counter, but cannot be shorter.)

Line 2 of the code doesn’t trust you to stop calling it when you’re already done. This is an interesting hedge, but not really an oath. It may rise to that level at some point though if enough code starts to depend on this utterance.

Lines 3 & 4 are the only uses of (a non-nil) counter. They work just fine for any sequence. As such, while this function always returns either nil or an array, it will accept any sequence as the counter.

Line 5 outright promises that the function is going to return nil. This is, of course, countermanded on line 7. I suppose this is akin to verses 3, 4, & 5:

[3] Or if a woman makes a vow to the LORD, and binds herself by some agreement while in her father’s house in her youth, [4] and her father hears her vow and the agreement by which she has bound herself, and her father holds his peace, then all her vows shall stand, and every agreement with which she has bound herself shall stand. [5] But if her father overrules her on the day that he hears, then none of her vows nor her agreements by which she has bound herself shall stand; and the LORD will release her, because her father overruled her.

Even though the loop on line five is released of its promise, it still feels like a cheap back door.

So, you can see that even in these eight lines of code, I have made many, many inadvertent promises. Hopefully, in future code, I will keep a more careful tongue.

TC Lispers, July Presentations Online July 19th, 2009
Patrick Stein

The July meeting of the Twin Cities Lisp Users Group was this past Tuesday. There were four presentations on the agenda:

The presentation slides and videos of the talks are available through the links above, and directly at the tclispers.org site. Enjoy!

Lisp Jump Tables July 12th, 2009
Patrick Stein

A few weeks back, I posted about how I know that I have more to grok about Lisp macros. In that post, I needed to explain Lisp macros a bit for those not in the know. The example that I used was creating a jump table.

Geoff Wozniak just posted an article in his blog wherein he explores using a macro to do jump tables using an array of functions. He performed timing tests using Tim Bradshaw’s ncase macro versus case and switch.

The timing differences there are quite amazing. I’ve toyed around with optimizations quite a bit and have never seen Allegro outperform SBCL before this. It is also interesting that CLISP may already do this jump table optimization. I have never seen CLISP come anywhere close to SBCL or Allegro before this.

Update: Including timing numbers from my own machine

My machine is a 2.4GHz Intel-based Mac with 4GB of RAM. Here are the numbers that I get for this test. (All numbers in milliseconds).

  n case ncase switch
SBCL 1.0.29 10 166 149 207
  100 185 146 191
  1000 572 147 207
CLISP 2.47 10 115 211 188
  100 117 205 219
  1000 111 206 189
ACL 8.1 (heap-limited) 10 10 20 90
  100 20 10 90
  1000 1420 20 130
CMUCL 19f 10 120 120 180
  100 150 120 210
  1000 520 120 210
Clozure v1.3-r11936 10 19 15 128
  100 108 14 130
  1000 1214 14 117

  • With Allegro (ACL), I couldn’t compile easily. I had to load it, then compile it, then load the compiled version. If I just started Allegro and tried to compile, it complained that WITH-GENSYMS was an undefined function.
  • With a heap-limited LispWorks, I failed to even compile the program. *shrug*

Updates In Email

Email:

RSS Feeds

l