fhtr

Shooting high ISO in broad daylight

2012-01-29T05:22:00.001+02:00

DSLRs these days get usable results even at super-high ISO (sensor sensitivity to light, higher it is the less light you need). But, um, what's the use of ISO 25600 if most of your shooting happens in bright daylight (or anything other than pitch-black darkness). Let's think!

What you get from cranking up ISO is the capability to use faster shutter speeds and tighter apertures. So you get photos with less motion blur and more of the image in focus. So you could shoot in motion and not need to focus. This is starting to sound useful. Shoot while walking without having to stop.

How about framing though? Having your camera up on your eye when you're walking up stairs sounds like a recipe for broken bones. Shooting from the hip would be nicer, but now you can't see what's in frame. Wide-angle lenses to the rescue! Get enough of the scene in the image that you'll likely have your subject in the frame, and do the framing in post.

It's really fast to take photos if you don't have worry about focus or framing. Point camera at what you're interested in, press shutter, done.

High ISO looks like crap in color though. Go black and white and you'll get smooth usable results at 6400 and noisy results at 25600 on a Nikon D5100 / D7000 / Sony NEX-5N (all have the same sensor AFAIK. I have a D5100). I'd kinda like to try a NEX-5N with a pancake lens for a small setup.

To recap: set ISO to 6400 or 25600, shoot in black and white, use manual focus (set to near-infinity or 20 meters or something), set shutter speed to 1/1000 or 1/500, aperture to f/16, use a 24mm lens or wider, snap away while walking!

Here's a gallery of my results from yesterday. They're not all pure examples of this technique, for some I brought the camera to my eye to do framing.

Animating a million letters with WebGL

2012-01-27T14:24:00.003+02:00

Here's an WebGL animated book demo! It's got just 150000 letters, but it does scale up to two million.

Writing efficient WebGL is a bit funny. The basic idea is to collect as many operations as possible into a single draw call, as changing the WebGL state machine state and doing WebGL calls is relatively expensive. If you want to draw more than a couple thousand objects at once, you need to adopt a quite different strategy for drawing.

The usual way of drawing with WebGL is to set up your uniforms, buffers and shaders for each object, followed by a call to draw the object. Unless your object is very complex, the time taken in this way of drawing is dominated by the state setup. To draw in a fast way, you can either do some buffer editing in JavaScript, followed by re-uploading the buffer and the draw call. If you need to go even faster, you can push more computation to the shaders.

My goal in this article is to draw a million animated letters on the screen at a smooth framerate. This task should be quite possible with modern GPUs. Each letter consists of two textured triangles, so we're only talking about two million triangles per frame.

Ok, let's start. First I'm going to create a texture with the letter bitmaps on it. I'm using the 2D canvas for this. The resulting texture has all the letters I want to draw. Then I'm going to create a buffer with texture coordinates to the letter sprite sheet. While this is an easy and straightforward method of setting up the letters, it’s a bit wasteful as it uses two floats per vertex for the texcoords. A shorter way would be to pack the letter index and corner index into one number and convert that back to texture coordinates in the vertex shader.

I also upload a two-million triangle array to the GPU. These vertices are used by the vertex shader to put the letters on the screen. The vertices are set to the letter positions in the text so that if you render the triangle array as-is, you get a basic layout rendering of the text.

With a simple vertex shader, I get a flat view of the text. Nothing fancy. Runs well, but if I want to animate it, I need to do the animation in Javascript. And JavaScript is kinda slow for animating the six million vertices involved, especially if you want to do it on every frame. Maybe there is there a faster way.

Why yes, we can do procedural animation. What that means is that we do all our position and rotation math in the vertex shader. Now I don't need to run any JavaScript to update the positions of the vertices. The vertex shader runs very fast and I get a smooth framerate even with a million triangles being individually animated every frame. To address the individual triangles, I round down the vertex coordinates so that all four points of a letter quad map to a single unique coordinate. Now I can use this coordinate to set the animation parameters for the letter in question.

The only problem now is that JavaScript doesn’t know about the particle positions. If you really need to know where your particles are, you could duplicate the vertex shader logic in JavaScript and update them in, say, a web worker every time you need the positions. That way your rendering thread doesn’t have to wait for the math and you can continue animating at a smooth frame rate.

For more controllable animation, we could use render-to-texture functionality to tween between the JavaScript-provided positions and the current positions. First we render the current positions to the texture, then tween from the JS array towards these positions, updating the texture on each frame. The nice thing about this is that we can update a small fraction of the JS positions per frame and still continue animating all the letters every frame. The vertex shader is tweening the positions.

Using MediaStream API

2012-01-27T14:14:00.002+02:00

Using the MediaStream API to access webcam from JavaScript:

navigator.webkitGetUserMedia("video,audio",
function(stream) {
var url = webkitURL.createObjectURL(stream);
videoTag.src = url;
videoTag.onerror = function() {
stream.stop();
alert('camera error');
};
},
function(error) {
alert(error.code);
}
);

Very basic math

2012-01-27T04:06:00.002+02:00

I was playing around with the idea of presenting fractions in the same way as negative numbers. Instead of 1/x, you'd write /x. Just like instead of 0-x, you write -x. And since multiplication with single-letter symbols is often annotated with putting the symbols next to each other, marking the inverse with /x looks quite natural: A x /B = A/B, 9 x /7 = 9/7.

It also makes you think of the inverse in less magical terms. Consider the addition rule for fractions:

A C AD BC AD + BC
- + - = -- + -- = -------
B D BD BD BD

There's some crazy magic happening right there. The literal meaning is (A x D x 1/B x 1/D) + (C x B x 1/D x 1/B), but you wouldn't know from looking at that formula. And it gets even more confusing when you start multiplying and dividing with fractions. Think about the following for a moment:

A C AD
- / - = --
B D BC

Right?

In linear notation with /B and /D and suchlike, this all actually sort of makes sense in a non-magical way. Here's the first of the above two examples (with intermediate phases written out):

(A x /B) + (C x /D)
= [1 x (A x /B)] + [1 x (C x /D)]
= [(D x /D) x (A x /B)] + [(B x /B) x (C x /D)]
= [(A x D) x (/B x /D)] + [(B x C) x (/B x /D)]
= (/B x /D) x [(A x D) + (B x C)]

[here's where you go: "oh right, /7 x /4 = /28", analogous to 7 x 4 = 28]

And the second one:

A x /B x /(C x /D)
= A x /B x /C x D
= (A x D) x (/B x /C)

Note the similarity with addition:

A + -B + -(C + -D)
= A + -B + -C + D
= (A + D) + (-B + -C)

Now, you might notice that there is a bit of magic there. How does /(C x /D) magically turn into (/C x D)? Or -(C + -D) to (-C + D) for that matter. Let's find out! Here's how it works:

/(C x /D)
= 1 x /(C x /D)
= [(/C x D) x /(/C x D)] x /(C x /D)
= (/C x D) x /(/C x C x D x /D)
= (/C x D) x /(1 x 1)
= (/C x D) x /1 -- Remember the axioms 1 x N = N and N x /N = 1. Since 1 x /1 = 1 we get /1 = 1.
= (/C x D) x 1 = (/C x D)

For the -(C + -D) case, replace / with -, x with + and 1 with 0.

And there you have it, my small thought experiment. And derivations for some basic arithmetic rules. I kinda like how breaking the magic bits down into the basic field axioms makes things clearer.

[edit]

Why is /A x /B = /(A x B)?

/(A x B) x (A x B) = 1

1 x (/A x /B) = (/A x /B)
/(A x B) x (A x B) x (/A x /B) = (/A x /B)
/(A x B) x (A x /A) x (B x /B) = (/A x /B)
/(A x B) x 1 x 1 = (/A x /B)
/(A x B) = (/A x /B)

Fast code

2012-01-21T18:55:00.000+02:00

I was thinking of the characteristics of high-performance language runtimes (read: execution times close to optimal for hardware) and came up with this list:

Flat data structures (err, like an array of structs where struct N is in memory right before struct N+1)

streaming memory reads are prefetcher-friendly, spend less time chasing pointers

Tightly-packed data

memory fetches happen in cache line -sized chunks, tightly-packed data gives you more payload per memory fetch
fit more payload into cache, faster subsequent memory accesses

Reusable memory

keep more of the working set of data in cache

Unboxed values

spend less time chasing pointers
generate tight code for data manipulation because data type known (float/double/int/short/byte/vector)

Vector instructions

more bytes manipulated per instruction, data moves faster through an execution unit

Parallel execution

more bytes manipulated per clock cycle, data moves faster through the processor

Keep data in registers when possible

less time spent waiting for caches

Keep data in cache when possible

less time spent waiting for memory
instead of going over the full data set several times end-to-end, split it into cache-sized chunks and process each one fully before moving onto the next one

Minimize amount of unnecessary data movement between processing units

keep data close to processor until you're done with it, even more important with GPUs

Flat code layout

low amount of jumps per byte processed

Tight code

keep more of the program in cache

Interleaved I/O

work on already loaded data while loading in new data
minimum amount of time spent waiting for I/O

You might notice that the major theme is optimizing memory use. I started thinking of program execution as a way to read in the input data set and write out the output data set. The sizes of the input and output data give you a nice optimum execution time by dividing the data set size by memory bandwidth (or I/O bandwidth if you're working on big things). The flow of the program then becomes pushing this river of data through the CPU.

Suck in the data in cache line -sized chunks, process entire cache line before moving to the next, preload next cache line while processing the current one. Use vector instructions to manipulate several bytes of data at the same time, use parallelism to manipulate several streams of data at the same time. Make your processing kernel fit into L1 instruction cache

Gnnnn ok, back to writing JavaScript.

Favicon notify

2011-11-08T18:16:00.001+02:00

Hack of the day from a few weeks back: show small notification bubble in the favicon. Check out the code at https://github.com/kig/faviconNotify and go to http://fhtr.org/faviconNotify for a demo.

Usage


  FaviconNotify.set(number);
  FaviconNotify.clear();

Longer example of use in HTML


  <html>
    <head>
      <link rel="icon" href="favicon.ico">
      <script src="faviconNotify.js"></script>
      <script>
        window.onblur = function() {
          FaviconNotify.set(1);
        };
        window.onfocus = function() {
          FaviconNotify.clear();
        };
      </script>
    </head>
  </html>

Riffing on the work of +Michael Mahemoff http://softwareas.com/dynamic-favicons
And http://userscripts.org/scripts/show/24430
And +Mathieu Henri http://www.p01.org/releases/DEFENDER_of_the_favicon/

See also http://faviconist.com/favicon-library (updated dynamic favicon library to incorporate badge-setting).

Updated Three.js deck

2011-11-08T17:54:00.002+02:00

I updated the "Basics of Three.js"-presentation to work with the latest version of the library. My Google Developer Day "Introduction to WebGL"-presentation is also online, and updated to match the latest changes to Three.js.

In more detail, the changes were:

$FOO.addChild and $FOO.addLight were merged into a single method, $FOO.add
Camera was deprecated and split into a PerspectiveCamera and OrthographicCamera
Camera no longer has a target that it tracks, you need to use camera.lookAt(vec3) on every frame to accomplish that.
MeshShaderMaterial was deprecated and renamed to ShaderMaterial
material.ambient now needs an AmbientLight in the scene to work (AmbientLight color is multiplied by material.ambient in the shader)
ColladaLoader was renamed to THREE.ColladaLoader

Basics of Three.js Presentation

2011-09-27T17:54:00.002+03:00

The slides for my Basics of Three.js talk are now up at fhtr.org/BasicsOfThreeJS

They take you through getting started with Three.js, building small apps and using custom shaders. There's a lot of live examples to play with, including the world's worst 3D modeler (pictured). You can fork the deck at github.com/kig/BasicsOfThreeJS

If you want to show off your awesome Three.js skills, add a new section to the deck and send me a pull request. Remember to add your name to the credits if you do!

Canvas image filter library

2011-08-29T02:20:00.002+03:00

Here's a small image filter library for the canvas element: github.com/kig/canvasfilters

It's based on the HTML5Rocks Canvas Image Filters article.

WebGL Filesystem Visualizer

2011-08-27T15:23:00.005+03:00

I put the HTML5WOW 3D filesystem visualizer demo online at fhtr.org/wfsv. Works only on Chrome as it uses the file input webkitdirectory attribute to recursively list directory contents.

Use the file input button in the top-left corner to select a directory to visualize. The visualizer tries to generate a model of the full directory tree, so it only works on small dir trees (less than 500 files or so). Use the console to navigate by clicking on file names. You can also use some shell commands like ls, cp, rm and mv. The filesystem contents aren't modified though, only the visualization.

Here's a video of the visualizer in action from the Google I/O 2011 "HTML5: The Wow and the How"-presentation:

Ken Burns effect using CSS

2011-08-19T15:54:00.006+03:00

The Ken Burns effect is a special effect used in documentaries when you only have a static photograph of an interesting item. To add some movement and life to the photograph, you zoom into the photo and pan towards a point of interest. It's named the Ken Burns effect because it was used a lot by a documentary film maker named Ken Burns.

Anyhow.

You can achieve the Ken Burns effect using CSS animations. It's not even particularly difficult. Just create a div with overflow:hidden to hold the image, then change the image's CSS transform property. Or if you want to be totally retro and backwards-compatible, you could also achieve the effect by changing the image's top, left, width and height using a JS setInterval.

So, CSS:


.kenburns {

  overflow: hidden;

  display: inline-block;

}

.kenburns img {

  transition-duration: 5s;

  transform: scale(1.0);

  transform-origin: 50% 50%;

}

.kenburns img:hover {

  transform: scale(1.2);

  transform-origin: 50% 0%; /* pan towards top of image */

}

And the corresponding HTML:


<div class="kenburns" style="width:640px; height:480px;">

  <img src="image.jpg" width="640" height="480">

</div>

If you hover over the image, it will slowly zoom in and pan towards its top edge.

You can see the effect in action here. And a more complex version with a JS-driven lightbox here.

The (quick, hacky) code is on GitHub.

Spinner library

2011-08-11T22:27:00.003+03:00

Ok, uploaded the spinner library to GitHub. There's a demo page too. It's a wee bit buggy though. I think the transition events are not firing when the page is hidden.

Even getting it to the current phase was a pain, the interaction between CSS transitions and animations, the DOM, JavaScript and events is flaky at the edges. Once I get an animation or transition going, it works well. But if I want to change some CSS properties in JS with the intent to trigger an animation or transition, it's pretty much guesswork whether it will actually work without glitches.

Hmm... As a workaround, I think I could just create a new element on every show/hide-cycle. That way the animation can't possibly fire one 100% frame before starting from 0%, right?? (This happened on Firefox... or maybe it was an unanimated frame or something. I ended up switching from animations to transitions because of that.)

On Chrome Canary, border-radius and overflow:hidden don't behave like you'd expect. Instead of clipping the contents to the rounded box, it clips on the rectangular box. I don't know why. And sometimes the rendering flickers when using transforms that toggle HW compositing for the page.

Story of Tomte's New Hat

2011-07-19T00:24:00.003+03:00

To break the monotonic programmingification, here's a story about Tomte's new hat.

Microbenchmark findings

2011-07-06T12:34:00.002+03:00

Looping through canvas ImageData pixels: 1D traversal for (var i=0; i<data.length; i+=4) ... is slightly faster than 2D traversal for (var y=0; y<height; y++) { for (var x=0; x<width; x++) ... }. Caching the ImageData width, height and data properties to local variables is slightly faster than not (except in the 2D case, where it's a good bit faster). Stay away from crazy counting hacks, they're more likely to slow you down than speed you up (probably because they make the compiler's job harder and it can't optimize the code properly).

Based on the above, my preferred pixel loop format is:


var width = id.width;
var height = id.height;
var data = id.data;
for (var y=0; y<height; y++) {
  for (var x=0; x<width; x++) {
    var off = (y*width + x) * 4;
    var r = data[off];
    var g = data[off+1];
    var b = data[off+2];
    var a = data[off+3];
    // ...
  }
}

And if you just need to map over the pixels and don't care about the coords, the simple 1D loop:


var data = id.data; // or just use id.data directly, no big perf impact
for (var i=0; i<data.length; i+=4) {
    var r = data[i];
    var g = data[i+1];
    var b = data[i+2];
    var a = data[i+3];
    // ...
}

Sparse blitting vs. dense blitting: If all you need to do is change the values of a small number of pixels, use fillRect. If you need to fill the entire canvas or a significant portion (say, 256x256) of it with custom data, use putImageData. There's no major performance difference between clearing an ImageData buffer before use vs. allocating a new one, so I'd go with clearing the buffer just to avoid extra GC work.

Text drawing: fillText is significantly faster than strokeText. Whether the text is aligned on an integer pixel only seems to matter on Opera's fillText.

Path drawing: Nothing very conclusive, path point count doesn't seem to have a major impact on path point throughput.

Spritesheets: For best performance, cache your spritesheet frames to separate canvases.

Image drawing: Align your images to integer pixels, don't transform them, use canvas elements instead of IMGs. Transformations perform very badly on non-accelerated canvases. Just offsetting an image by fractional pixels causes the browser to use the slow path. On accelerated canvases, transformations don't really matter all that much. You still get the best perf by doing aligned non-transformed draws though.

Clearing the canvas: Just use clearRect.

WebGL texture sources: Use premultiplied textures if possible. Use ImageData if possible, just don't specify it as premultiplied (same goes for typed arrays). Canvases are faster than images on Chrome, images are faster than canvases on Firefox. Typed arrays are about as fast as canvases.

New articles

2011-06-08T16:53:00.002+03:00

I have two new articles up at HTML5 Rocks. The first is about making image filters using the canvas element and the second one is about building a kiosk-style app for showcasing Chrome Experiments at this year's Google I/O.

Go check them out and let me know how you like them.

Canvas & WebGL microbenchmarks

2011-06-02T16:09:00.003+03:00

I've been writing some microbenchmarks for canvas and WebGL over the past few weeks. They're at jsPerf, but jsPerf being what it is, I'm going to collect all the benchmark links here. So, without further ado:

Canvas 2D tests

Clearing the canvas
Drawing images
Spritesheets vs. individual sprites
Path creation
Filling and stroking text
Drawing text at different sizes
Sparse blitting
Dense blitting
Looping through ImageData pixels

There's also this test on the speed of the different canvas composite operations over at GitHub, I probably should port it over to jsPerf as well.

WebGL tests

Texture sources and pixelStorei settings

Runfield & Remixed Reality

2011-04-24T00:10:00.010+03:00

Oh right, I did these demos for Mozilla a couple months back:

Runfield

Runfield is a Canabalt clone with painted graphics (I made a guide on how to do your own graphics for it, but it's a bit "First you sketch a good-looking picture and then you paint it! Done!"). The graphics were painted in MyPaint and GIMP. The renderer is done with Canvas 2D and uses drawImage to draw thin vertical slices from the background images to make up the undulating ground.

The main things I wanted to communicate with Runfield were speed and polish. Showing that you can make a fast 2D game with JS and have it look good. Accordingly, most of the dev time was spent making the graphics and optimizing the engine (:

Optimization tips: draw images aligned to the pixel grid, eliminate overdraw (if you know that a part of an image is not going to show, don't draw that part), use the first couple seconds to detect the framerate and drop down to a lighter version if the framerate is low.

Remixing Reality

Remixing Reality is another demo to showcase what you can do with JavaScript today. It's processing video frames in real-time to locate AR markers and uses WebGL to draw 3D models on top of the markers. If you click the play button on the side, music starts playing and there's a 3D music visualizer powered by BeatDetektor2 and the Mozilla audio API, again analyzing the audio in real-time.

The AR library powering the thing is JSARToolKit, a pure-JS port of the Flash FLARToolKit (which uses NyARToolKitAS3, which is a port of the Java NyARToolKit, which is a port of the C ARToolKit. Whew.) Porting it over to JS was pretty quick, since the AS3 syntax is close enough to JS syntax that I could write a good-enough syntax translation script in a couple days. Then I implemented the AS3 class semantics in JS and off we go.

Well, it wasn't quite that easy. The syntax translator is a hack and I had to go and manually fix things. And implement the pertinent parts of Flash's BitmapData. And write a shim to make it work with Canvas. But hey, 14 kloc port in a week!

The job didn't end there though. It was slow. The biggest slowdown was that the library was reading a couple pixels at a time from the canvas, and each of those reads called getImageData. So, cache it, problem solved.

It was still a bit slow, mostly due to FLARToolKit using BitmapData's color bbox queries to do feature detection. I.e. find the smallest rectangle in the bitmap that includes all pixels of a certain color. Each call to BitmapData#getColorBoundsRect needs to go through the pixels in the bitmap and find the first row where the wanted color is found, then the bottom row, then scan the rows in between to find the left-most and right-most columns. This process was not too fast in JS.

But NyARToolKit, the library which FLARToolKit is based on, was doing the feature detection in an entirely different way. Its algo was running on RLE-compressed images (Run-Length Encoding: pack data as [value][number of repetitions], e.g. aaabbbb becomes a3b4). And since the images in question are thresholded to black and white, RLE works very well. Expected result: smaller images => less work for JS => faster.

So... I made the JS version use the NyARToolKit version. And hey, it was 5x faster! Nice!

Another thing that helped performance on Firefox 4 was using typed arrays instead of normal JS arrays. Fx4's JIT generates more efficient machine code for typed arrays. On Chrome 10(? IIRC), typed arrays and normal arrays didn't have much of a performance difference, but the code ran fast enough on normal arrays already.

For the 3D stuff I used my Magi library. With a Blender export script to get the models in. And a slightly tweaked lighting shader to make it fill the unlit areas with some ambient. Fun times.

Browser rendering loop

2011-04-01T02:19:00.003+03:00

The browser rendering loop is how the browser displays the web page to you.

The main stages of the rendering loop are:

Updating the DOM.
Rendering the individual elements.
Compositing the rendered elements to the browser window.
Displaying the browser window to the user.

The DOM updates happen in JavaScript or in CSS transitions and animations. They include things like "Hey, I'd like that header to turn red." and "Please draw a thick line on the canvas."

If you're drawing to a canvas, you might expect the browser to draw as soon as you issue a drawing command. Which is what actually happened in earlier browser versions. But in the latest browsers, it doesn't quite work that way. Nowadays the browser queues up the drawing commands and only starts drawing when it absolutely needs to. Which is usually just before compositing, in the second stage of the rendering loop.

However, if you want to force the browser to finish drawing before continuing JS execution, you can try doing getImageData on the 2D Canvas and readPixels in WebGL. As they need to return the finished image, they should force the browser to flush its draw queue. This comes in handy if you ever need to figure out the time it took for the browser to execute your drawing commands.

Once all the individual elements are drawn, the browser composites them together to create the final browser window image. And finally, the browser window image is shown to the user via the OS window manager.

The frame rate perceived by the user is the frequency at which step 4 is repeated. In other words, how often the updated browser window is shown to the user.

As most flat panel displays can only update 60 times per second (the TV frequency), the browser tries to display only up to 60 frames per second. Going over 60 when your display can't take advantage of it would only burn more CPU and reduce battery life, so it makes sense to clamp the update frequency to the display's update frequency.

Optimally, the browser would finish all its drawing before doing a new composite, but current browser implementations have slight problems with that. So it's really rather difficult to figure out the actual framerate visible to the user. If you have a high-speed video camera, you could record the display and see how fast it's updating.

But if you want to do it all in the browser, you could try something like this. First, hook up to the frame loop with requestAnimationFrame. Second, flush the drawing queue when your frame is done. Third, measure time from flush to flush. Hopefully browsers will move towards requestAnimationFrame only being called after flushing the previous frame.

Firefox and Chrome actually make this a bit easier for you by providing some built-in framerate instrumentation. Chrome dev channel has an about:flags FPS counter that (sadly) only works for accelerated content. Firefox 4 has the window.mozPaintCount property that keeps track of how many times the browser window has been redrawn.

References:
GPU Accelerated Compositing in Chrome
Hardware Acceleration in the latest Firefox 4 beta
ROC: Measuring FPS
Measuring HTML5 Browser FPS, or, You're Not Measuring What You Think You're Measuring

Please send me a note if anything above is wrong / misguided / an affront to your values and I'll try and fix it.

Detecting new globals in JavaScript

2011-03-31T21:37:00.003+03:00

JavaScript has this annoying little feature where if you forget the 'var' keyword when assigning to a variable, it creates a new global. This has a tendency to cause some veeery interesting bugs. So. It would be nice to detect those new globals. And I happen to have just the thing for that.

Paste this snippet to your webpage (preferably after creating all the globals you intended to create):

<script type="text/javascript">
  if (true /* MONITOR_GLOBALS */) {
    (function(){
      var globals = {};
      var startGlobals = [];
      for (var j in window) {
        globals[j] = true;
        startGlobals.push(j);
      }
      if (false /* PRINT_INITIAL_GLOBALS */)
        console.log("Initial globals: "+startGlobals.sort().join(', '));
      setInterval(function() {
        var newGlobals = [];
        for (var j in window) {
          if (!globals[j]) {
            globals[j] = true;
            newGlobals.push(j);
          }
        }
        if (newGlobals.length > 0)
          console.log("NEW GLOBALS: "+newGlobals.sort().join(', '));
      }, 1000);
    })();
  }
</script>

Now whenever a script creates a new global, you get a notification in your JS console. Hopefully sparing you from some agonizing hours of debugging.

Trying to make sense of this whole cooling thing

2011-02-05T20:22:00.007+02:00

I kind of find CPU heat sinks interesting. Probably because of the noisy fans messing with my concentration. Why do we need those things anyhow?

I'm not a physicist so if I'm saying Very Silly Things below, I'd appreciate if you left a comment to set me straight.

Fans are devices to move ambient-temperature coolant fluid (read: air) into contact with heat sink fins, so as to maximize heat conduction from the heat sink to the environment. Heat conduction = thermal conductivity * area * temperature difference / distance between the different temperatures ^{[according to Wikipedia]}. If you think of the heat sink surface–ambient air -system in terms of that model, a fan would be trying to minimize the distance between the hot heat sink and the cool air. In practice, pushing off the boundary layer of hot air and replacing it with cooler air.

There are also other methods for optimizing the heat conduction equation, here's a small roundup. Endothermic reactions for soaking up heat and transporting it away: heat pipes. Maximize the temperature difference: dry ice, LN2 (plus the boiling is endothermic). Increase thermal conductivity of the coolant fluid: water cooling, mineral oil cooling. Increase conducting area: finned heat sinks, larger heat sinks - and I guess the increased mass acts as a thermal buffer that soaks up spikes in heat production. Decrease distance between the different temperatures: fan pushing in cold air, chimney enhancing natural convection. (Natural convection: hot air expands, making it less dense than cool air, which lets cool air fall below hot air. As long as you have gravity, that is.)

Sci-fi, where our Author reveals his lack of understanding in The Art of Physicks

You could also get cold air close to the heat sink fins by ionizing air and charging the heat sink with reverse polarity to attract the cool air to the heat sink (hopefully neutralizing the air in the process). Maybe that could get cold air more effectively through the boundary layer than mechanical pushing with a fan? They're doing some commercial research into that, and there's also a DIY Ion Cooler.

Some sort of channels on the heat sink surface to require less airflow to replace the boundary layer. Funnel the air to a higher velocity to push through the boundary layer more efficiently.

This is neat: vibrating piezoelectric heat sink fins.

Heat travels as phonons, maybe you could somehow create a heat waveform and sap it out through destructive interference with an audio source. Sort of like noise cancelling headphones. Audible noise is in the kilohertz range, heat is in the gigahertz range. Apparently these guys at MIT are making phononic mirrors.

Make the heat do work, cooling the heat sink down in the process. Stuff a heat sink full of piezoelectric crystals to convert the thermal expansion of the heat sink into electricity, conduct the electricity away to cool the heat sink.

Can heat be focused? Heat up a small region of the heat sink to a high enough temperature for blackbody radiation to really kick in, use mirrors to transport the resulting photons away. (Intensity of blackbody radiation grows as fourth power of absolute temperature, temperature grows roughly linearly with power input (apart from phase changes).

A meter-high 0.1x0.1m chimney with 50C internal air and 20C room temperature could generate 0.01 m^3/s airflow, or 22 CFM? Plugging in values to Q = C*A*sqrt(2*g*H*(Ti/T0-1)): 0.7 * 0.01m^2 * sqrt(2*9.81m/s^2*1m*(323K/293K-1)) = 0.0099 m^3/s.

Rotate the chimney, get a fire tornado!

An interesting page on thermal design

More half-baked ideas to keep the year rolling

2011-02-05T11:24:00.002+02:00

Tone down flashing ads by rendering the web page to an accumulation buffer at a low opacity, so that there'll be a motion blur effect that evens out the flashing.

Intel X25-V SSD

2011-01-27T10:34:00.006+02:00

Bought an Intel X25-V 40GB SSD to do some random read benchmarking. Don't stick it into a PATA controller driven SATA port, as that limits the bandwidth to 70 MB/s and random access times to 5 ms.

Sequential read speed on the unused disk was limited by the SATA bus. So I guess it's not actually reading anything for unused blocks, just generating a string of zeroes on the controller and sending it back. Would be pretty amusing to have a special disk driver that keeps a list of the unused blocks in system RAM and generates the data on the CPU for any read accesses to them. You could probably keep it all in L1 and get nice bogus 4k random read benchmark numbers. "1 ns access latency!?? 100 GB/s random read bandwidth??? Whaaaat?" (An allocation bitfield for 80 gigs in 512k blocks is about 20 kB. Increase blocksize or add RLE for bigger volumes.)

Read speed on actual data was around 200 MB/s. Average random read time for a 4k block was 0.038 ms with 128 threads and 0.31 ms with a single thread. So the controller can do 8 reads in parallel. Which is nice for a cheapo SSD.

Didn't really test write performance apart than doing a cursory check. And yes, it is low bandwidth. 40 MB/s streaming writes. So it's best used as a random read drive.

Might make a nice random read array with 6 drives or somesuch. Hypothetically: 48 parallel reads, 160 000 4k reads per second (650 MB/s). And then all you need is software that can take advantage of that.

Getting flash chip latencies lower would be good.

Stupid ideas to kick off 2011

2011-01-17T04:00:00.004+02:00

They have those spinning disks with a couple gigs of flash as read cache, right. So, how about spending about 30e to put a bunch of fast flash chips and a couple gigs of DRAM on a motherboard to act as SATA read cache. Read the flash into DRAM during boot checks (I guess you have around 5 seconds to do it, so 0.5-1GB/s should suffice.) Ooh, mysteriously computer boots as if from ramdisk.

Or alternatively, if you like software more than hardware, use the fast flash as the OS disk, suck it to RAM during boot checks, write OS to be able to utilize that.

Also, maybe they could make a computer that's not actually a computer but a rock. And then you could have a mouse that's not a mouse but a rat and it'd bite your fingers off and give you rabies.

How to take advantage of that memory?

2011-01-15T10:38:00.012+02:00

[edit] Here's an rc.d init script I'm using to prefetch the FS to page cache: usr-prefetch. Copy to /etc/init.d and run `update-rc.d usr-prefetch start 99 2 .` to add it to runlevel 2 init scripts. I currently have 4GB of RAM, and the total size of the prefetched stuff is 3.7GB, so it might actually help. Maybe I should go buy an extra 2GB stick. It takes about a minute to do the prefetch run with cold cache, so the average read speed is around 60 MB/s. Which is pretty crappy compared to the 350 MB/s streaming read speed, maybe there's some way to make the prefetch a streaming read.

Ubuntu's ureadahead reads in the files accessed during boot in disk order. If you add all your OS files to it, it should be possible to stream the whole root filesystem to page cache in about ten seconds. I dunno.

And now I'm reading through ureadahead's sources to figure out how it sorts the files in disk order. It's using fiemap ioctl to get the physical extents for the file inodes, then sorts the files according to the first physical block used by the file. To read the files to the page cache, it uses readahead. If there was some way to readahead physical blocks instead of files, you could collapse the physical blocks into spans with > x percent occupancy, do a single streaming read per span (and cache only the blocks belonging to the files to cache).

Another easy way to do fast caching would be to put all system files on a separate partition, then cache the entire partition with a streaming read. Or allocate a ramdisk, dd the root fs onto it at boot, use it through UnionFS with writes redirected to the non-volatile fs. But whether that's worth the bother is another thing.
[/edit]

In the previous post I put together a hypothetical machine with 16GB RAM and 1.5GB/s random access IO bandwidth. Which is total overkill for applications that were built for machines with 100x less resources. What is one to do? Is this expensive hardware a complete waste of money?

The first idea I had on taking advantage of that amount of RAM was to preload the entire OS disk to page cache. If the OS+apps are around 10GB in total, it'll take less than 10s to load them up at 1.5GB/s. And you'll still have 6GB unused RAM left, so the page cache isn't going to get pushed out too easily. With 10s boot time, you'd effectively get a 20GB/s file system 20 seconds after pushing the power button.

Once you start tooling around with actual data (such as video files), that is going to push OS files out of page cache, which may not be nice. But what can you do? Set some sort of caching policy? Only cache libs and exes, use the 1.5GB/s slow path for data files? Dunno.

If you add spinning disks to the memory hierarchy, you could do something ZFS-like and have a two-level cache hierarchy for that there filesystem. 10GB L1 in system RAM, 240GB L2 on SSDs, a couple TB on spinning disks. But again, how sensible is it to cache large media files (which will likely be the primary use of the spinning disks.) IT IS A MYSTERY AARGH

It'd also be nice to have a 1-2GB of faster RAM between the CPU caches and the system RAM (with GPU-style 150GB/s bandwidth), but maybe you can only have such things if you have a large number of cores and don't care about latency.

But hmm, the human tolerance for latency is around 0.2s for discrete stuff, 0.015-0.03s for continuous movement. In 0.2s you can do 4GB worth of memory reads and 0.3GB of IO, in 0.02s the numbers are 0.4GB and 0.03GB respectively. On a 3.5 GHz quad-core CPU, you can execute around 3 billion instructions in 0.2s, which'd give you a ratio of 0.75 ops per byte for scalar stuff and 6 ops per byte for 8 op vector instructions. On GPUs the ops/byte ratio is larger, at close to 20 ops/byte for vector ops, which makes them require more computation-heavy algorithms than CPUs.

In 4GB you can fit 60 thousand uncompressed 128x128 images with some metadata. Or hmm, 6 thousand if each has ten animation frames.

The computer hardware of 2011

2011-01-13T23:00:00.007+02:00

Here's a 1500e computer for 2011. It's a build where I tried to minimize the bottlenecks and get a nice memory pyramid with GPU cache 2 TB/s, GPU RAM 300GB/s, CPU cache 200GB/s, RAM 20GB/s, IO 2GB/s. Throw in 100e worth of spinning disks there in RAID-0 for a trailing 0.2GB/s disk IO.

CPU: 4-6 cores, 3.5GHz, 150e.
Motherboard: 150e.
Case: 50e steel box with razor-sharp edges and PSU from hell.

Total: 350e.

Memory: 8GB @ 150e.
Memory bandwidth: 20GB/s.

Total: 500e.

IO subsystem: 4xSSD RAID-0 @ 500e.
2GB/s streaming IO bandwidth.
1.5GB/s random access IO bandwidth at 4k block size.

Total: 1000e.

GPU: two upper middle-class things @ 500e.
Computing power: 3 TFLOPS.
Memory bandwidth: 300GB/s.

Total: 1500e.

(Optionally 2x 500GB disks, 0.2GB/s IO bandwidth @ 100e for a total of 1600e.)

Reasons for the selections:

The performance difference between a 150e CPU and a 300e CPU is ~15%, and around 25% with overclocking. The point is more to get high memory bandwidth with minimal cost. Non-bargain mobo selected for 6Gbps SATA and dual GPUs. Only 8GB memory to cut costs. All money saved was tossed on getting random access IO onto the order-of-magnitude curve and the GPUs to cap the top end of the memory bandwidth curve (and for a disproportionate amount of computing power). The SSD numbers are based on recently announced SATA 6Gbps drives. If you don't care about GPU performance, get a 300e CPU, use IGP, double the RAM and buy an extra SSD.

Now, programming this thing is going to be interesting. Not only do the CPU and GPU require parallelism to extract reasonable performance, but the random access IO subsystem does too. The IO subsystem latencies are probably in the 0.2ms range but it can do 64 parallel accesses.

A single-threaded C program will achieve maybe 1% of peak performance of the hardware, if that.

For legacy software, you could buy a 300e computer and it'd perform just as well.

CPU: 70e, mobo: 60e, case: 50e, RAM: 20e, SSD: 90e. Toss a 40e disk in there for good measure.

Or better yet, quit computers altogether and spend all the money on furniture.

The Space Campaign

2010-12-26T14:42:00.004+02:00

Space Captain, First Galactic Baron, Royal Ring Knight, First Solar King and Third Moon Prince of the Space Planet New New Zealand, Second Master of the Holy Star Order, Imperial Great Lord of Space and Independent Space General, Sergeant Major Lieutenant Stanislaw Stanisgreat stood up from his space chair on the space deck to stare across the spacebow of the space ultrabattledestroyer HMS Stiff Upper Lip with his spacescope. There! In the distant depths of space, near the space planet of New South Macedonia, a blacker than black space cutter was cutting space towards the space merchant fleet sailing the space waves towards the space port of New North America. The weak, defenseless merchantmen were sitting in their space deckchairs just waiting to be picked off from space by the evil boat of space.

Space Captain Stanislaw cursed a terrible space curse under his space breath to no one in particular, which was just as well for there is no sound in space. He gave the order to turn the space wheel of the space ultrabattledestroyer towards the space cutter, hoist the main space sail to the space mast and blow the space horn to sound the space crew to readiness. The space sails billowed in the space wind and started the ultrabattledestroyer towards the nefarious space cutter at infinite velocity. Even at infinite velocity, Stanislaw knew he'd be too late to stop the space cutter from sinking the space merchantmen into the depths of space and feeding their space crews to the space fishes.

An endless instant later, the space ultrabattledestroyer closed to engagement range from the space cutter and turned its right side towards the dastardly space cutter. Stanislaw ordered a full thirty-fold space broadside alpha strike to be fired at the space cutter at his mark. Mark! The crewmen of the space cutter shouted in great alarm at the ultrabattledestroyer but their shouts were never heard in space before the report of the space guns of the ultrabattledestroyer drowned them out under a terrible booming fury of space.

Millions of space cannonballs flew across space in slow motion, tearing gaping space holes in the space sails of the space cutter. One space cannonball hit the cutter right in the middle and split it in two. The mangled space cutter turned its cowardly tail and tried to run across space to the space asteroid field in the middle of space to escape the space pursuit of the space ultrabattledestroyer.

Too bad for the space cutter, the ultrabattledestroyer had its space sensors deployed in space, giving them full knowledge of space. Space Captain Stanislaw set his magnificent space plan in motion and deployed a single space rowboat right in the space path of the space cutter. The space rowboat was carrying a deadly cargo of space gunpowder, which the spacemen aboard it set to explode with a slow-burning space fuse. The spacemen jumped to space off the space rowboat and space swam back to the welcoming space arms of the space ultrabattledestroyer.

The space captain of the space cutter, Shining Diamond Vizier Democratically Elected Big Cat Corporate Chief Communist Politruck Liberal Nazi Warlord Space President of Badmanistan, General Colonel of Three Stripes and One Star Badstar Badmanis noticed the space rowboat a mere picosecond too late. His bellowed order to space duck was cut short in space as the space rowboat exploded with thunderous space roar and shot three million nuclear space lightnings everywhere in space, shattering the space cutter into five trillion pieces. An endless shower of space blood colored the space red. That was one space cutter that would not cut space again, thought the entire space crew of the space ultrabattledestroyer in unison in space.

Space Captain Stanislaw Stanisgreat turned to face his space crew and congratulated them for a space job well done. His lily-white space suit and space tophat had been splattered crimson from the space blood and his visage was like that of a butchered space goat. The Space Admiral Peerless Thirty-fold Queen Emperor, Defender of Space, Attacker of the Realm, First Dame of Spacistan, Royala Bigshotikins climbed the ladder aboard the ultrabattledestroyer from her million-mile long gigaspacehyperdreadnoughtcarrier space flag ship and bestowed the space medal of Biggest Imperial Glory of Space upon the beaming Space Captain Stanislaw Stanisgreat. Stanislaw fought back a single space tear of infinite boundless space pride and saluted for hours at the Space Admiral.

THE END
THESE BOOKS MAKE ME ANGRY PUBLISHING, LTD.

WebGL Orrery

2010-12-19T08:56:00.004+02:00

I wrote a small WebGL orrery a month back or so, using planet and orbit data from Wikipedia. The orbits came out a bit wrong (planet positions on orbits are wrong, the planet rotation axis rotates along the orbit, planets velocities are constant, the orbits / axes may well be flipped somewhere), but it's a passably "realistic" orrery. Well, realistic apart from the whole having to scale down orbits by a factor of thousand to make the planets visible -thing. With 1:1 orbits, the planets are smaller than a pixel in size when viewed from other planets and the sun is about a pixel in diameter when viewed from Neptune.

It's quite extensible too, it should be possible to add in all the moons just by filling in their data. The coordinate system and z-buffer aren't accurate enough to do anything fancy like a spaceship sim with 1:1 orbits though :)

A deck of cards

2010-12-01T14:57:00.012+02:00

If not for the order property, you could cram a deck of cards into a 64-bit int. Though if you had a RNG cyclical in full_deck_size, that'd work to preserve the order property. Which would be enough for poker and other stackless or write-only stack card games.

With the cyclical RNG you could pop a card off a full deck in O(1). With incomplete decks (that is, a deck with non-RNG-generated cards removed) you'd need to go through random cards until the generated card is in the deck. The worst-case time-complexity for popping all cards off an incomplete deck would be O(full_deck_size) and the average complexity for popping one card off an incomplete deck would be O(full_deck_size/incomplete_deck_size). Anyhow, incomplete decks would be as fast or slow as full decks, with an extra existence check slowing them down a bit further. The RNG would need a seed with log₂(fac(full_deck_size)) bits, so log₂(52!) = 226 bits for a deck without jokers.

A naive version of the RNG would be to generate a random number between 0 and 52!-1, then convert the number into 52 digits in factorial base. Divide first by 51!, the result is the first card. Divide the remainder by 50!, the result is the second card (note that you have to remove the first card from consideration for the second card index). Continue dividing the remainders until you have 52 cards. Or just generate a shuffled 52 element array of numbers and go through that :P

... I should get back to work now.

Source


#include <inttypes.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>

typedef int64_t card_mask;
typedef int8_t card;

typedef struct deck {
  card cards[55];
  int8_t length;
  card_mask mask;
} deck_t;

#define suit_size 13

static const card_mask hearts = suit_size*0;
static const card_mask spades = suit_size*1;
static const card_mask clubs = suit_size*2;
static const card_mask diamonds = suit_size*3;
static const card_mask jokers = suit_size*4;

static const char* suits[5] = {"Hearts", "Spades", "Clubs", "Diamonds", "Jokers"};
static const char* ranks[13] = {
  "Ace", "Two", "Three", "Four", "Five", "Six", "Seven", "Eight", "Nine", "Ten",
  "Jack", "Queen", "King"
};

card_mask remove_cards_from_mask(card_mask deck, card_mask cards)
{
  return deck & ~cards;
}

card_mask add_cards_to_mask(card_mask deck, card_mask cards)
{
  return deck | cards;
}

card_mask query_mask(card_mask deck, card_mask cards)
{
  return deck & cards;
}

card_mask mask_of_card(card c)
{
  return (card_mask)1 << (card_mask)c;
}

card card_for_suit_and_rank(card_mask suit, card_mask rank)
{
  return suit + rank;
}

int deck_initialize(deck_t *deck, int jokers)
{
  card c = 0;
  if (jokers < 0 || jokers > 3) 
    return 0;
  deck->length = 52 + jokers;
  deck->mask = ((card_mask)1 << (card_mask)deck->length)-1;
  for (c=0; c<deck->length; c++) {
    deck->cards[c] = c;
  }
  return 1;
}

void deck_clear(deck_t *deck)
{
  deck->length = 0;
  deck->mask = 0;
}

void deck_shuffle(deck_t *deck)
{
  int8_t i,j;
  card tmp;
  for (i=deck->length-1; i>0; i--) {
    j = rand() % (i+1);
    tmp = deck->cards[j];
    deck->cards[j] = deck->cards[i];
    deck->cards[i] = tmp;
  }
}

card_mask deck_query_mask(deck_t *deck, card_mask m)
{
  return query_mask(deck->mask, m);
}

card_mask deck_query_card(deck_t *deck, card c)
{
  return deck_query_mask(deck, mask_of_card(c));
}

card_mask deck_query_deck(deck_t *deck, deck_t *deck2)
{
  return deck_query_mask(deck, deck2->mask);
}

int deck_add(deck_t *deck, card c)
{
  card_mask m = mask_of_card(c);
  if (query_mask(deck->mask, m)) {
    return 0;
  } else {
    deck->mask = add_cards_to_mask(deck->mask, m);
    deck->cards[deck->length] = c;
    deck->length++;
    return 1;
  }
}

int deck_remove(deck_t *deck, card c)
{
  int8_t i=0, j=0;
  card_mask m = mask_of_card(c);
  if (!query_mask(deck->mask, m)) {
    return 0;
  } else {
    deck->mask = remove_cards_from_mask(deck->mask, m);
    for (i=0,j=0; i<deck->length; i++) {
      if (deck->cards[i] != c) {
        deck->cards[j] = deck->cards[i];
        j++;
      }
    }
    deck->length = j;
    return 1;
  }
}

card deck_pop(deck_t *deck) 
{
  card c;
  if (deck->length < 1) {
    return -1;
  } else {
    deck->length--;
    c = deck->cards[deck->length];
    deck->mask = remove_cards_from_mask(deck->mask, mask_of_card(c));
    return c;
  }
}

const char* suit_string(card s)
{
  if (s >= suit_size*5) 
    return "Invalid";
  return suits[s / suit_size];
}

const char* rank_string(card s)
{
  return ranks[s % suit_size];
}

void print_card(card c)
{
  printf("%s of %s\n", rank_string(c), suit_string(c));
}

void print_deck(deck_t *deck)
{
  int8_t i;
  printf("Deck %p:\n", (void*)deck);
  for (i=0; i<deck->length; i++) {
    print_card(deck->cards[i]);
  }
}

void print_mask(card_mask m)
{
  int i;
  const char* symbols[5] = {"♥", "♠", "♣", "♦", "☻"};
  for (i=0; i<64; i++) {
    if (i%suit_size == 0) printf(" %s", symbols[i/suit_size]);
    if ((m>>i) & 1) putchar('#');
    else putchar('.');
  }
  putchar('\n');
}

int main() 
{
  int i=0;
  deck_t deck, hand;
  card c;
  srand(time(NULL));
  deck_initialize(&deck,2);
  deck_initialize(&hand,0);
  deck_clear(&hand);
  deck_shuffle(&deck);
  printf("Deck %p has %d cards and mask %ld.\n", (void*)&deck, deck.length, deck.mask);
  print_mask(deck.mask);
  for (c=0; c<54; c++) {
    if (!deck_query_card(&deck, c)) {
      printf("Error: Deck doesn't have a card.\n");
      print_card(c);
    }
  }
  printf("Hand %p has %d cards and mask %ld.\n", (void*)&hand, hand.length, hand.mask);
  print_mask(hand.mask);
  for (c=0; c<54; c++) {
    if (deck_query_card(&hand, c)) {
      printf("Error: Hand has a card it shouldn't have.\n");
      print_card(c);
    }
  }
  deck_remove(&deck, hearts);
  if (deck_query_card(&hand, hearts)) 
    { printf("Error: Hand has a card it shouldn't have.\n"); print_card(hearts); }
  deck_remove(&deck, spades);
  if (deck_query_card(&hand, spades)) 
    { printf("Error: Hand has a card it shouldn't have.\n"); print_card(spades); }
  deck_remove(&deck, clubs);
  if (deck_query_card(&hand, clubs)) 
    { printf("Error: Hand has a card it shouldn't have.\n"); print_card(clubs); }
  deck_remove(&deck, diamonds);
  if (deck_query_card(&hand, diamonds)) 
    { printf("Error: Hand has a card it shouldn't have.\n"); print_card(diamonds); }
  printf("Deck %p has %d cards and mask %ld.\n", (void*)&deck, deck.length, deck.mask);
  print_mask(deck.mask);
  for (i=0; i<5; i++) {
    c = deck_pop(&deck);
    print_card(c);
    deck_add(&hand, c);
  }
  printf("Deck %p has %d cards and mask %ld.\n", (void*)&deck, deck.length, deck.mask);
  print_mask(deck.mask);
  printf("Hand %p has %d cards and mask %ld.\n", (void*)&hand, hand.length, hand.mask);
  print_mask(hand.mask);
  print_deck(&hand);
  if (deck_query_deck(&deck, &hand) || deck_query_deck(&hand, &deck)) {
    printf("Error: Deck and Hand intersect.\n");
  } else {
    printf("Deck and Hand are disjoint.\n");
  }
  for (i=0; i<hand.length; i++) {
    if (deck_query_card(&deck, hand.cards[i])) {
      printf("ERROR ");
    } else {
      printf("OK ");
    }
  }
  printf("\n");
  for (i=0; i<deck.length; i++) {
    if (deck_query_card(&hand, deck.cards[i])) {
      printf("Error: Found a card in Deck that is in Hand.\n");
      print_card(deck.cards[i]);
    }
  }
  print_card(suit_size*5);
  return 0;
}

Output


$ ./a.out 
Deck 0x7fff98907250 has 54 cards and mask 18014398509481983.
 ♥############# ♠############# ♣############# ♦############# ☻##..........
Hand 0x7fff98907210 has 0 cards and mask 0.
 ♥............. ♠............. ♣............. ♦............. ☻............
Deck 0x7fff98907250 has 50 cards and mask 18013848686551038.
 ♥.############ ♠.############ ♣.############ ♦.############ ☻##..........
Eight of Hearts
Seven of Clubs
Ten of Diamonds
Six of Diamonds
Nine of Hearts
Deck 0x7fff98907250 has 45 cards and mask 17714777228828286.
 ♥.######..#### ♠.############ ♣.#####.###### ♦.####.###.### ☻##..........
Hand 0x7fff98907210 has 5 cards and mask 299071457722752.
 ♥.......##.... ♠............. ♣......#...... ♦.....#...#... ☻............
Deck 0x7fff98907210:
Eight of Hearts
Seven of Clubs
Ten of Diamonds
Six of Diamonds
Nine of Hearts
Deck and Hand are disjoint.
OK OK OK OK OK 
Ace of Invalid

Yet more stupid ideas

2010-11-23T17:35:00.004+02:00

Amateur astronomer telescopes, live video streaming, combine the different observations into high-res frames: live 24/7 HD webcams of the planets.

Speckle imaging
Compressed sensing
Image Deblurring with Blurred/Noisy Image Pairs

More stupid ideas

2010-11-23T11:52:00.003+02:00

Take The European Citizens' Initiative, add a public forum, link it up with EU proceedings, build a Citizens' Parliament.

The ECI proposal says that you need a million signatories to send a citizens' initiative to the European Commission. The EC is the main proposal-generating entity in the EU. The European Parliament, in contrast, is the proposal-filtering entity (now with new proposal-generating powers granted to it by the Lisbon Treaty). In that sense, the ECI impacts the very top of the EU legislative system.

The million signatories must include 0.2% of the voting-age citizens for a third of the member states to be valid. And after collecting 300 000 statements of support, you must submit it to the EC for an admissibility acid check where the EC figures out whether they can legally do it or not. After two months, you should get back a yay/nay. And you can collect the statements of support online as long as your mechanism is secure (further details pending).

After that, 3 months of checking the signatures for validity by the member states using random sampling. After that, 4 months of examination by the EC, at the end of which it "would then be required to set out its conclusions on the initiative and the action it intends to take in a communication, which would be notified to the organiser as well as to the European Parliament and Council and would be made public."

Not the most blisteringly fast method of democratic participation, is it? Still, it's very interesting.

I wonder how much governance is a problem of thinking power (in terms of "need more heads to think better solutions!") and how much a sensor network problem ("need to know what's wrong with the world to fix it!") and how much a problem of integrating the two into a working system. How would this kind of a wide-ranging massively large parliament work best?

A large and diverse population has experts from all the different walks of life. Which is probably where the thinking power of the citizens' parliament would best manifest. But you need to find the best proposals from the noise of mediocre proposals. Then you have the sensor network side, where a large population is again an asset. But again, there's a large amount of noise with a large population. So perhaps the important thing is making a good network to filter the signal.

Stupid ideas

2010-11-18T10:06:00.005+02:00

#1
Use toString, Hungarian notation and a JS parser to compile JS functions to GLSL. To run the JS shaders in JS, implement the GLSL built-ins as a JS library and write a pure-JS GL implementation to laugh at.


// Pseudo-code time!

var myShader = {
  vertex: {
    attributes: {v3Vertex: buffer, v2TexCoord: anotherBuffer},
    uniforms: {m4PMatrix: cameraMatrix, m4MVMatrix: worldMatrix},
    varyings: ['v2TexCoord0'],

    // the hypothetical JS GLSL runtime would do something like 
    // shader.void_main.apply(GLSLContext, [])
    void_main: function() {
      var v4v = this.vec4(this.v3Vertex, 1.0);
      this.v2TexCoord0 = this.v2TexCoord;
      this.gl_Position = this.m4PMatrix.mulV4(this.m4MVMatrix.mulV4(v4v));
    }
  },

  fragment: {
    varyings: ['v2TexCoord0'],
    uniforms: {s2DTexture0: myTexture, fRadius: 0.1},
    void_main: function() {
      var fd = this.fRadius / 5.0;
      var v4Color = this.vec4(0.0);
      for (var fi=0.0; fi<5.0; fi++) {
        var v2tc = this.v2TexCoord0.add(this.vec2(fi*fd,fi*fd));
        v4Color = v4Color.add(this.texture2D(this.s2DTexture0, v2tc));
      }
      this.gl_FragColor = v4Color;
    }
  }
};

puts(GLSL.convert(myShader.vertex));
=>
precision highp float;
attribute vec3 v3Vertex;
attribute vec2 v2TexCoord;
uniform mat4 m4PMatrix;
uniform mat4 m4MVMatrix;
varying vec2 v2TexCoord0;
void main() {
  vec4 v4v = vec4(v3Vertex, 1.0);
  v2TexCoord0 = v2TexCoord;
  gl_Position = m4PMatrix * (m4MVMatrix * v4v);
}

puts(GLSL.convert(myShader.fragment));
=>
precision highp float;
varying vec2 v2TexCoord0;
uniform sampler2D s2DTexture0;
uniform float fRadius;
void main() {
  float fd = fRadius / 5.0;
  vec4 v4Color = vec4(0.0);
  for (float fi=0.0; fi<5.0; fi++) {
    vec2 v2tc = v2TexCoord0 + vec2(fi*fd, fi*fd);
    v4Color = v4Color + texture2D(s2DTexture0, v2tc);
  }
  gl_FragColor = v4Color;
}

Maybe you could do an OpenCL implementation like that too.

#2
3D software should come with sculpt-friendly pre-rigged model kits with UV maps (yes, having some trouble with the technical aspects of 3D modeling). Also, voxel sculpt looks amazing.

Storage visions

2010-10-22T06:45:00.002+03:00

Some sort of user-friendly UI to ZFS+Plan9. One big main filesystem that spans several drives and automatically does striping, mirroring and versioning. Fast drives are used as cache, big drives as mirrors and history.

Plug in a new internal drive and usable space, redundancy and total performance grows by some percentage of the new drive's specs. Bring a new computer in the local area network and the filesystems on the computers are merged into a network-wide filesystem, adding usable space, redundancy and performance.

Pull out a drive and the file system rebalances itself over the remaining drives. If some files require the pulled out drive, display notification and freeze the files until the drive is restored. Use a disk utility to separate a drive from the filesystem into a new filesystem, the utility copies files unique to the drive to other drives and disconnects and formats the separated drive when done.

Computers on the local network are used as page cache and as mirrors of the filesystem. Cloud service is used as a backup mirror and remote sync. Treat the entire memory hierarchy and network as parts of the same single filesystem. When on the network, it doesn't matter which computer or which disk has the file, apart from performance.

Transient devices carry a small local cache of the important things on the filesystem, and access the rest of the filesystem over the network. Either by connecting to the home network or by going through the cloud service. The local cache of the filesystem prioritizes metadata over the actual data and access to the actual data can usually be done by streaming a low-bandwidth version.

*mirrorshades*
*techno soundtrack*
*sunrise*

How many dwarfs live in Mt. Everest

2010-10-20T03:55:00.006+03:00

Dwarfs, those earth-dwelling little buggers, grand masters of volume. What if they burrowed the entire Mount Everest into a big dwarf-warren? Let's find out!

First, some assumptions.

Dwarfs are clearly good at cramped living as they live in mostly closed-off subterranean cities. They have also mastered the art of stonecraft to a magical level, so their buildings may well compare favorably to modern ferro-concrete towers. Let's assume that a dwarven city boasts a similar population density as urban Hong Kong, i.e. 250000 dwarfs per square kilometer with an average building height of 45 meters. A cubic kilometer of dwarvopolis would consist of 20 such layers and have a population of 5 million dwarfs.

Mount Everest is kinda pointy, but not all that pointy. Let's suppose that its slope is 45 degrees and that it approximates a cone. That would give it around 720 cubic kilometers in volume. But don't forget that dwarfs are subterranean, so we could reasonably assume that the dwarven city of Sagarmatha extends an extra kilometer under the surface with dwarven suburbia sprawling in a 10 kilometer radius. With the subterranean suburbs added, our dwarven city would have a total volume of around one thousand cubic kilometers.

From which a simple multiplication leads us to a total Sagarmatha dwarvopolitan population of 5 billion.

But why stop at Everest? Cast your net wider! The entire population of the Great Himalayan Dwarven Realm would be... the area of Himalayas is around 610000km^2, Tibetan plateau height is 4.5km, assume 0.5% of Everest population density for the sticks, i.e. 25000 dwarfs per cubic kilometer, and put it all together for a grand total of 69 billion Himalayan dwarfs.

Ubuntu 10.10 nits

2010-10-20T00:23:00.008+03:00

Nice fonts.
Nice boot screen.
Nice boot screen doesn't work with fglrx.
Dark theme bg & window style screws up labels with dark text.
Dark title bar & dark menu bar occasionally look silly with the light-grey window bg. Especially Firefox looks bad. The thing with dark bg themes is that the amount of apps designed to work properly with one is damn close to zero. (Hence, I suppose, the light-grey window bg.)
The active window titlebar gradient looks kinda sloppy. Ditto for popup menu gradients.
Flat window buttons walk a tightrope between sucking and being just boring.
GTK widget theme buttons are very nice.
Dropdown widget has the triangle vertically off-center by 1px?
The salmon active color in the widget theme is not very nice.
Scrollbar widget looks different from the rest of the theme and a bit weird. Exacerbated by rounded buttons but no alpha channel, so Firefox textarea scrollbar buttons look icky.
New Rhythmbox icon is nice.
Rhythmbox itself is still nowhere near Amarok 1.4.
Amarok 2.0 is still nowhere near Rhythmbox.
Gnome Terminal icon is vertically off-center in top panel.
Notification area in top panel looks nice with the monochrome theme.
App icons in top panel are too tall and look bad.
Weather applet rainy icon has a black border. Other widget icons (other weather icons too) don't have borders, making the rainy icon look out-of-place.
The speech bubble icon in top panel is one pixel too low (more like one pixel too high).
The networking applet icon is super-confusing. I keep looking at it and wondering "Wtf was that icon for again? Check updates? No.. oh right, networking."
The (very low-visibility) round corners in panel dropdown menus look nice. But suffer from lack of text padding.
Window title color is actually orange instead of dark grey (try Chromium with GTK theme).
1px wide window resize handles make me cry.
The purple color works surprisingly well.
Default desktop bg is nice.
Popup text labels have too little horizontal padding. Otherwise they look nice.
Date & time applet text has too little horizontal padding as well.
The panel applets have crap spacing so you need to manually move them around by a couple pixels to make the layout look less horrible.
PulseAudio is glitchy and the volume control does not work (you can only set it between mute, medium volume and super-loud volume, no low volume setting available).
The new fancy volume control doubling as a music player remote is nice.
udev screws up md RAID arrays by being clueless (starts only one disk of a RAID-1 -> stops boot process thx to unable mount disk & dirty array after reassembly -> need resync -> FFFFUUUUUUU. How about using mdadm -A --scan like sane people do?) [edit: mdadm -A --scan wants to have an md device name like /dev/md0 instead of /dev/md/MyRAID, I dunno if that sunk udev or if it's just doing silly things.]
Nautilus is slow so I very much try not to manage files with it (instead terminal for quick & programmable things, Filezoo for overviews and visual sorting)
The CPU speed indicator applet could show all the CPUs at once instead of requiring running several copies of it.
I don't get what the Social Networks thing does. It doesn't seem to do anything. Maybe that's what it does.
What is Ubuntu One? Perhaps these things should be explained as one starts them up?
Why is this Gwibber thing so slow? Or perhaps it just lacks a progress indicator, making it feel slow.
Evolution feels like a really fiddly and slow interface to Gmail.

I do really like it for the most part, at least after a year in Debian's permabroken Gnome-land. I'd use FVWM or something cool and edgy and configurable like that for window management but compiz no worky on non-retarded WMs and theming FVWM (while infinitely easier than theming metacity) is limited and painful. Not that I've done any theming for a couple years. Plus you don't have a nice GUI to make F1 start a new terminal (just a nice thousand-line config file).

Maybe there should be a WebKit-based desktop environment.

Fast SSD is fast

2010-10-18T01:00:00.007+03:00

The PCI-E SSD I bought is a 120GB OCZ RevoDrive. Which is basically two 60 GB OCZ Vertex 2 SSDs mounted on a PCI-E fakeRAID controller.

Works nicely in Ubuntu 10.10, though apparently less nicely on earlier versions (problems booting from fakeRAID).

The amusing thing about the SSD is that the controller is compressing all the data it writes to the drive. So you get 540MB/s speeds for reading a file full of zeroes but just 300MB/s for reading a file full of MP3. Write speeds are 475MB/s for zeroes, 300MB/s for MP3s.

A random access 4kiB read takes around 0.26 ms, for a 16MB/s random read speed. But you can do around 8 random reads in parallel, which gets you a 117MB/s random read speed. My 7200 RPM disks can do around 75MB/s with a streaming read, so 117MB/s random access speed is absolutely nuts.

Parallel reads are kinda hard to program though. You could spawn a bunch of threads and do a gather into a shared array I suppose, though that feels a bit heavy (hah, right. I expect the threading overhead to be negligible.) Do the reads from a shared mmap of the file, each thread reading its own segments off it? [edit: no, mmap apparently serializes its page faults or OpenMP doesn't run the threads in parallel]

Little benchmark

2010-10-16T17:46:00.003+03:00

Got my computer put together yesterday, and now's the time to benchmark it! On a Saturday afternoon...

Image correlation algorithm benchmark:


OpenCL
  240 GBps -- Athlon II X4 640, 3GHz (12GHz aggregate), 2MB L2
   85 GBps -- Core 2 Duo E6400, 2.1GHz (4.3GHz aggregate), 2MB L2 
OpenMP+SSE optimized
  103 GBps -- Athlon II X4
   45 GBps -- Core 2 Duo
OpenMP+SSE naive
   13 GBps -- Athlon II X4
    5 GBps -- Core 2 Duo

Pretty much linear scaling with clock frequency in OpenCL. Both have a 3 cycle L1 latency and the algorithm is very much an L1 cache benchmark, so this isn't too surprising. The SSE version has some bandwidth / load-balancing bottleneck going on, and the naive version is pretty much a pure memory bandwidth benchmark.

New hardware!

2010-10-13T00:49:00.002+03:00

In our cheap (and not-so-cheap) upgrade category, got me an Athlon II X4 + mobo + RAM (300€) and a PCI-E SSD (300€). Which should get me roughly 3x more CPU power, double the memory bandwidth, 5x the IO bandwidth and something like 700x random access performance. Now just have to piece it together and benchmark!

WebGL presentation editor prototype

2010-10-10T13:20:00.004+03:00

Try it out here, requires a WebGL capable browser. It's a bit buggy, so watch your head.

You can check out the sources on GitHub. What's still missing is ~~making the editor work properly without WebGL~~ DONE (though it's not Mobile Safari -friendly), making the slide editing more efficient, and ~~fixing some markup parser bugs~~ fixed some. At that point it'd be a decent barebones slideshow app. Add themes, fancy filters and animations and stir on low heat.

Desk drawer coding

2010-10-04T15:17:00.006+03:00

Tally of different sorts of crap I've done to pass the time:

Music players written: 4 (with mplayer / <audio> / Flash for playback, so maybe these are more "playlist players").
Image slideshows written: 5, three of which on arbitrarily long lists.
Video players /w playlist: 1 (I just repurposed one slideshow).
Document readers: 3, though they just show images of pages.
Scene graph libraries for different retarded purposes: 7.
GUI toolkits: 1. But it didn't go too well.
LaTeX-style text justifier: 1.
Layout text inside a polygon: 1. And it seems like there's zero literature on text layout algorithms.
Tilemap zoomers: 3.
Infinite tree hierarchy zoomers: 1.
File managers: 1. And 3 web-based file managers.
File thumbnailer libraries: 2.
Thumbnail cache systems: 3.
OpenGL-based presentation apps: 3.
File annotation systems: 1.
Metadata editing forms: 3.
File metadata extraction libraries: 1.
Standard library rewrites: 1.
Array combinator libraries: 2. Well, maybe 3.
Search query parsers: 2, the second time Doing It Right with a parser generator.
Search query matchers: 1.
DOM helper libraries: 2.
Context-sensitive help in REPL: 1.
Distributed message-passing libraries: 2.
Network service discovery libraries: 1.
Media players with networked storage and playback (controller on mobile, storage across the network, playback on HTPC): 1, but durr.
Converter graph with pathfinder to convert files: 1 (IIRC plugged to the media player too, so it played PDFs with text-to-speech).
Animation libraries with tweens and timelines: 1.
SQL databases (tables and indices, not database engines): 5? Not very much.
Memcached to cache results of expensive queries: 1.
HTML fragment caches: 1.
Tar parsers: 1.
3D model parsers: 3 of varying levels of completeness.
Huffman encoder & decoder: 1 (yay, coursework).
Image file formats: One vector drawing hack, plus three that basically concatenate JPEGs into a bigger file.
Toy drawing programs: 3.
Games: 1.

Editors: 0.
Compilers: 0.
Desk calculators: 0.
Programming languages: 0 (unless you count an OO Brainfuck implementation).
Something useful: 0.

Hmm, that was more than I initially thought :<

Visitor stats for September 2010

2010-09-27T10:31:00.003+03:00

Browser market share

43% Firefox
31% Chrome
10% Safari
10% Explorer
4% Opera

OS market share

56% Windows
22% Linux
19% OS X
1.6% iOS
0.7% Android

More than 90% of visitors had a screen width of 1280 or higher.

94% had Flash, 1% had a Flash version other than 10. 85% had Java.

Most search engine visitors wanted to know about SSE matrix multiplication. Clearly a hot topic.

Composing shaders

2010-09-25T14:29:00.002+03:00

I haven't found a good way to do modular GLSL shaders, so here's some blather on that subject.

Shaders themselves are strings of source code, compiled to a binary by the OpenGL drivers. A shader has a void main() that is called to execute the shader. There are also global definitions in the form of attributes, varyings and uniforms. And a precision qualifier like precision highp float;, and C preprocessor macros.

A shader also has a type, it can either be a vertex shader or a fragment shader (or a geometry shader). The difference between the two is that a vertex shader determines what areas to consider for modification, whereas a fragment shader determines the output value for each element in the modified areas. Imagine you're drawing a filled circle: the vertex shader would output the bounding box of the circle and the fragment shader would test each element in the bbox to see if it lies inside the circle, setting them to the fill color if they do.

But back to the composing part. If you have some common helper functions that you'd like to use in your different shaders (e.g. defaultTransform(), blend functions), how do you reuse them? I tried to link different shader objects together but it didn't work for whatever the reason. So now I'm just concatenating strings.


ShaderSource = Klass({
  initialize: function(type, source) {
    this.type = type;
    this.text = toArray(arguments).slice(1).map(function(s) {
      return s.text ? s.text : s;
    }).join("\n");
  }
});

Maybe you want to further split your shaders into functional segments and turn those segments on and off for different uses. One way is to put some #ifdefs in the shader and prepend a block of #defines to your source.


ShaderSource = Klass({
  initialize: function(type, source) {
    this.type = type;
    this.text = this.origText = toArray(arguments).slice(1).map(function(s) {
      return s.text ? s.text : s;
    }).join("\n");
  },

  setDefines : function(defines) {
    this.defines = defines;
    var defs = [];
    for (var i in defines) {
      defs.push("#define "+i+" "+defines[i]);
    }
    defs.push(this.origText);
    this.text = defs.join("\n");
  }
});

Another way would be to affix an identifier to the source strings and only use a string if it's enabled. This is like writing a feature-poor preprocessor by yourself.


ShaderSource = Klass({
  initialize: function(type, source) {
    this.type = type;
    this.segments = toArray(arguments).slice(1);
    this.text = this.segments.map(function(s) {
      return s.text ? s.text : s;
    }).join("\n");
  },

  setEnabled : function(enabled) {
    this.enabled = enabled;
    this.text = this.segments.map(function(s) {
      if (!s.id || enabled[s.id])
        return s.text ? s.text : s;
      else
        return "";
    }).join("\n");  
  },

  setDisabled : function(disabled) {
    this.disabled = disabled;
    this.text = this.segments.map(function(s) {
      if (!s.id || !disabled[s.id])
        return s.text ? s.text : s;
      else
        return "";
    }).join("\n");  
  }
});

Or some mix of the above. I don't really know what's the super-best way to do this.

Easy inheritance in JavaScript

2010-09-25T06:28:00.005+03:00

JavaScript inheritance is cumbersome. You have a constructor function and it has a prototype. And inheritance is done by setting the prototype to a newly created parent object? What? That means that you better not do anything in the constructor, otherwise you'll have a royal pain deriving new prototypes from it.

So I did what I do and wrote my own inheritance system, based on mixins. It still functions pretty much like the JavaScript system, but is more convenient for building the constructor and the prototype from existing classes.

Here's the infamous OO geometry demo:


Rectangle = Klass({
  initialize : function(w, h) {
    this.w = w;
    this.h = h;
  },

  getArea : function() {
    return this.w * this.h;
  }
})


Square = Klass(Rectangle, {
  initialize : function(s) {
    Rectangle.initialize.call(this, s, s);
  }
})


Moving = Klass({
  mass : 1,

  initialize : function(n) {
    this.s = new Array(n);
    this.a = new Array(n);
    this.v = new Array(n);
    this.f = new Array(n);
    for (var i=0; i<n; i++){ 
      this.s[i] = this.a[i] = this.v[i] = this.f[i] = 0;
    }
  },

  move : function(t) {
    for (var i=0; i<this.s.length; i++) {
      var a = this.a[i], s = this.s[i], v = this.v[i], f = this.f[i];
      var newA = f/this.mass;
      var newV = v + a*t + (1/2)*(newA-a)*t*t;
      var newS = s + v*t + (1/2)*a*t*t + (1/3)*(newA-a)*t*t*t;
      this.a[i] = newA;
      this.v[i] = newV;
      this.s[i] = newS;
    }
    this.previousMass = this.mass;
  }
})


Printable = {
  valueString : function() {
    var a = [];
    for (var i in this)
      if (typeof this[i] != 'function') a.push(i);
    a.sort();
    for (var i=0; i<a.length; i++)
      a[i] = a[i] + ": " + this[a[i]];
    return a.join(", ");
  },

  print : function() {
    console.log(this.valueString());
  }
}


MovingRectangle = Klass(Rectangle, Moving, Printable, {
  density : 0.19,
  height : 1,
  initialize : function(w,h,d) {
    Rectangle.initialize.call(this,w,h);
    Moving.initialize.call(this,2);
    if (d)
      this.density = d;
    this.mass = this.getVolume() * this.density;
  },

  getVolume : function() {
    return this.getArea() * this.height;
  }
})

Klass is a function that builds a constructor. The returned constructor calls its prototype's initialize method.

The constructor's prototype is built by merging together all the arguments given in the Klass call. If an argument has a prototype property defined, the prototype is used for merging. Otherwise the argument is merged as itself.

The built prototype is then finally merged into the constructor function, so that you have access to the prototype functions and properties from the constructor object itself. Which is useful for doing calls to other classes, as you can do e.g. Parent.someMethod.call(this) instead of Parent.prototype.someMethod.call(this).


Klass = function() {
  var c = function() {
    this.initialize.apply(this, arguments);
  }
  c.ancestors = toArray(arguments);
  c.prototype = {};
  for(var i = 0; i<arguments.length; i++) {
    var a = arguments[i];
    if (a.prototype) {
      Object.extend(c.prototype, a.prototype);
    } else {
      Object.extend(c.prototype, a);
    }
  }
  Object.extend(c, c.prototype);
  return c;
}

Object.extend merges the src object's attributes with the dst object, ignoring errors. Ignoring the errors is useful in merging native objects (e.g. styles) as they've got some stuff that can't be overridden.


Object.extend = function(dst, src) {
  for (var i in src) {
    try{ dst[i] = src[i]; } catch(e) {}
  }
  return dst;
}

toArray creates a new array from an object that has a length property. Useful for dealing with all the silly Array-like objects such as the function arguments object and the lists returned by the DOM API.


toArray = function(obj) {
  var a = new Array(obj.length);
  for (var i=0; i<obj.length; i++)
    a[i] = obj[i];
  return a;
}

WebGL slides translated

2010-09-22T08:00:00.003+03:00

I translated my WebGL slides to English and put them on GitHub Pages. Which is great, btw. Upgraded my account to Micro just to give them props for that.

If you would like to buy an app to make fancy WebGL slideshows, send me a mail and I'll get back to you.

If you'd like to subscribe to a web-based zoomable file management & sync service instead, send me a mail and I'll get back to you.

In other news, I'm back on Twitter.

Saturday Night Links

2010-09-17T19:57:00.005+03:00

ARTtech 2010: Fairlight's rendering secrets from AssemblyTV. The render graph approach is interesting. And the demo editor.

Yay! SVG in IMG coming in Firefox. Should enable fancy scalable UI images and easy SVG textures in WebGL. (Edit: In 2010-09-18 nightly SVG in IMG works but can't be drawn to canvas or uploaded as WebGL texture.)

Creature workflow in Mudbox 2011.

Auto-retopology in 3D Coat. Take a high-poly model, paint over areas that you want at high resolution, draw lines to suggest edgeloop direction, done!

Maybe one day 3D modeling won't require thinking about polygons. You'd just have a high-resolution form that you deform and then save a low-res version for interactive use. I dunno. The whole "create box, extrude, cut, unwrap, texture"-flow turns me off. It's like you're doing reverse origami: you start with a paper box, then you extrude it, bend it, cut bits off, glue bits back on, cut it open and splay it flat, draw a picture on it, put it back together, stuff a wire model inside, animate. And all the time you're worrying if you got the bends just right so that the joints don't look stupid.

Frontend Finland WebGL slides

2010-09-16T22:22:00.009+03:00

I gave a talk at the Frontend Finland meetup about WebGL. Had a great time, thanks to the organizers! My presentation slides are available here in glorious Finnish. Best experienced with a browser that has a working WebGL implementation (fear not, there's a fallback to plain unadorned HTML). Writing the slideshow app (and the slides) took about three days, with four days of preliminary library work (and the app dev fed into the library too).

Here are some links on the topics covered:

The git repo for the utility library and the presentation can be found here.

The SVG presentation by Matti Paksula was really nifty as well. After the show, we were brainstorming about some kind of a super slideshow where you'd write the slides in HTML (with a nice CSS), then they'd get rendered using SVG / Canvas / WebGL depending on the capabilities of your browser. (What my presentation is doing is parsing the DOM into scene graph nodes, then rendering them using WebGL. Doing a similar system that rendered using SVG wouldn't be too difficult.)

Preserving the clickability of links and other such webness would be problematic with WebGL, but a mix of SVG, CSS 3D and WebGL might get you the best of both worlds: overlay SVG elems on top of the WebGL canvas, change the CSS 3D transform of the elems to match the WebGL scene graph. That'd be a lot like CAKE's mixed HTML/Canvas, except in 3D. If you want fancy FX on your SVG elems while preserving interactivity, you'd render them to textures and put the original element at opacity 0.01 on top of the fancy render, maybe? Or do shaders for SVG?

Saturday morning link

2010-09-11T04:11:00.002+03:00

Dust by Elude. Amiga AGA. Insane.

Dude!

2010-09-10T06:12:00.002+03:00

SDO video of solar activity

2010-09-09T19:37:00.002+03:00

SDO video of solar activity over almost three months (via VizWorld)

Arithmetic

2010-09-08T17:29:00.003+03:00

You have things, some things. Call these one-things. There is also a no-thing, call it 0. When you add one-thing to 0, we call that 1. Add a one-thing to 1 and you get 2. Add one-thing to 2 and you get 3. In that way we build a small list of numbers: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9.

When you add one-thing to 9, we don't have a number for that here. But let's call it a ten-thing and mark it on left of the one-things. We have 1 ten-thing and 0 one-things, written in short as 10. Then we can add a one-thing to the 0 one-things to get 11. When we add a one-thing to 19, we get two ten-things: 20. When we have ten ten-things, we call that a hundred-thing and write it as 100: one hundred-thing, zero ten-things, zero one-things.

In a similar fashion, we could have a shorter list of numbers, say, 0, 1, 2, 3, 4. Now when you add one-thing to 4, you get one five-thing and zero one-things, which we can write as 10₅. When we have five five-things, we call that a twenty-five-thing and write it as 100₅.

There are also negative things. When a negative thing is added to a positive thing, they both disappear. Adding one negative thing to one positive thing makes both disappear, leaving you with zero things. If you add 5 negative things to 8 positive things, the 5 negative things and 5 positive things disappear, leaving 3 positive things. Adding 2 negative things to 3 more negative things gets you 5 negative things. Adding 4 negative things to 2 positive things makes the 2 positive things and 2 negative things disappear, leaving 2 negative things.

When you add 5 negative one-things to 2 positive ten-things, you need to crack open one of the ten-things to get ten positive one-things. The 5 negative things and 5 positive one-things disappear, leaving you with 1 positive ten-thing and 5 positive one-things: 15.

To add 5 negative one-things to 1 positive hundred-thing, you first crack open the hundred-thing to get ten ten-things, then crack open one of the ten-things to get ten one-things. The 5 negative one-things and 5 positive one-things disappear, leaving you with 9 positive ten-things and 5 positive one-things: 95.

Adding a zero to a number results in the number itself. Adding 5 to 0 results in 5. Adding 0 to 5 results in 5.

Negating a number turns a positive number into a negative number and a negative number into a positive number. Positive 2 negated is negative 2. Negative 2 negated is positive 2. Adding a number to its negation is zero. Positive 2 added to negative 2 is zero.

We write positive numbers with a + in front of them and negative numbers with a - in front of them. Positive 2 is +2. Negative 2 is -2. We write addition with a + between the numbers. The + symbol is called plus. Negative 3 plus positive 2 is -3 + +2. We write negate as -(number). Positive 2 negated is -(+2). Negative 2 negated is -(-2).

Two numbers are equal if they are the same. Equality is written with a =. 2 equals 2 is written 2 = 2 and amounts to true. If the two numbers are not the same, the equality operation returns false. If both sides of the equality operation are identical, the equality operation returns true. 2 = 2 is true. 2 = 3 is false. 1 + 1 = 2 is undecided, rewriting it as 2 = 2 or 1 + 1 = 1 + 1 makes it true.

Subtraction is an operation where the first number is added to the negation of the second number. Subtraction is written with a - between the numbers. The - symbol is called minus. -3 minus +2 is written -3 - +2 and means -3 + -(+2).

Multiplication is an operation where each of the things in the first number is replaced by the second number. Multiplying +3 by -2 replaces each of the 3 things with a -2, giving you +(-2) + +(-2) + +(-2) = -6. Multiplying +15 by +3 requires splitting open the ten-things to replace each of the one-things inside them with a +3, giving you +15 +3-things, which amount to +45 when added together. Multiplying -2 by -5 replaces the two things with -5-things, giving you -(-5) + -(-5), which simplifies to +5 + +5 and amounts to +10. Multiplying +0 by -2 gives you zero -2-things, which amounts to zero. Multiplying -2 by +0 replaces the two things with +0, giving you -(+0) + -(+0) = -0 + -0 = -0.

Fractional things are parts of one-things, much like how one-things are parts of ten-things. Fractional things are defined by how many parts are needed to make up a one-thing. When you need two parts for one-thing, the fractional thing is called a half. When you need three parts, the fractional thing is one-third. With five parts needed, it's a one-fifth. Five one-fifths added together amount to one. Similarly three one-thirds amount to one. Fractional numbers are written as 1/parts. One-fifth is written +1/5. One-third is written as +1/3. Negative one-sixth is written -1/6.

Adding 1/6 to 1/6 results in 2/6. Adding 1/6 to 5/6 results in 6/6, which equals 1. Adding 1/6 to 1 is 1/6 + 6/6 = 7/6. Adding 1/5 to 4/5 results in 5/5, which also equals 1. Adding 1/5 to 1 is 1/5 + 5/5 = 6/5.

Multiplication is written with an x between the numbers. The x is called times. 3 times 5 is written 3 x 5 and means 3 multiplied by 5, amounting to 15.

Multiplying a number by one results in the number itself. 5 x 1 is 5. 1 x 5 is 5.

The inverse of a number multiplied with the number itself results in one. The inverse of a number is the fraction with number parts. Zero does not have an inverse. The inverse of +5 is +1/5. The inverse of -3 is -1/3.

Three times the inverse of three is 3 x 1/3 = 1/3 + 1/3 + 1/3 = 3/3 = 1.

Negative two times the inverse of negative two is -2 x 1/-2 = -(-1/2) + -(-1/2) = 1/2 + 1/2 = 2/2 = 1.

Division is an operation where you multiply the first number with the inverse of the second number. 3 divided by 5 is 3 x 1/5 and amounts to 3/5. 8 divided by 2 is 8 x 1/2 and amounts to 8/2, which adds up to 4.

Division is written with a / between the numbers. 5 divided by -3 is written 5 / -3 and means 5 x -1/3, and amounts to -5/3, or -1 + -2/3.

The allowed moves when manipulating an equation with additions are:

Apply associativity rule to two additions: a + (b + c) = (a + b) + c
Apply commutativity rule to an addition: a + b = b + a
Add the same value to both sides of the equality: a = a <=> a + b = a + b

For an equation with multiplications, the above moves apply as well, but there is an additional move when dealing with a mix of additions and multiplications:

Apply distributivity rule: a x (b + c) = a x b + a x c

These are the only allowed structural moves for manipulating equations of additions and multiplications.

Additions can use the following equalities for simplifying expressions: a + 0 = a and a + -a = 0. Multiplications can use a x 1 = a and a x 1/a = 1.

We can build more moves from the above, for example:


I'm playing fast and loose with the associativity in these, 
a + b + c = (a + b) + c = a + (b + c), 
I'm using whichever is the easiest at any given moment.

Also, multiplication binds tighter than addition:
a x b + c x d = (a x b) + (c x d)

     x + b = c      | a = a <=> a + b = a + b
x + b + -b = c + -b | a + -a = 0
     x + 0 = c + -b | a + 0 = a
         x = c + -b | a - b = a + -b
         x = c - b
x + b = c <=> x = c - b

      a x b = c       | a = a <=> a x b = a x b
a x b x 1/b = c x 1/b | a x 1/a = 1 when a != 0
      a x 1 = c x 1/b | a x 1 = a
          a = c x 1/b | a / b = a x 1/b
          a = c / b
a x b = c <=> a = c / b when a != 0

(a + b) x c = (a + b) x c   | a x b = b x a
(a + b) x c = c x (a + b)   | a x (b + c) = a x b + a x c
(a + b) x c = c x a + c x b | a x b = b x a
(a + b) x c = a x c + b x c

a x 1 = a x 1 | a x b = b x a
a x 1 = 1 x a

a x b = a x b                  | (a + b) x c = c x a + c x b
a x b = (a + -1) x b + (1 x b) | 1 x a = a
a x b = (a + -1) x b + b       | a - b = a + -b
a x b = (a - 1) x b + b

         -(a) = -(a)     | a = b <=> a + b = a + b
     a + -(a) = a + -(a) | a + -(a) = 0
     a + -(a) = 0        | a = b <=> a + b = a + b
-a + a + -(a) = -a + 0   | a + -a = 0
     0 + -(a) = -a + 0   | a + b = b + a and a + 0 = a
         -(a) = -a 

-1 x a = -1 x a     | 0 + a = a
-1 x a = 0 + -1 x a | a - b = a + -b
-1 x a = 0 - 1 x a  | 1 x a = a
-1 x a = 0 - a      | a - b = a + -b
-1 x a = 0 + -a     | 0 + a = a
-1 x a = -a

-1 x a = -a   | -a = -(a)
-1 x a = -(a)

-a x b = -a x b     | -1 x a = -a
-a x b = -1 x a x b | a x b = b x a
-a x b = a x -1 x b | -1 x a = -a
-a x b = a x -b

-a x b = -1 x a x b | -1 x a = -(a)
-a x b = -(a x b)

a x (b - c) = a x (b - c)        | a - b = a + -b
a x (b - c) = a x (b + -c)       | a x (b + c) = a x b + a x c
a x (b - c) = (a x b) + (a x -c) | a x b = b x a
a x (b - c) = (a x b) + (-c x a) | -a x b = -(a x b)
a x (b - c) = (a x b) + -(c x a) | a x b = b x a
a x (b - c) = (a x b) + -(a x c) | a - b = a + -b
a x (b - c) = (a x b) - (a x c)

a x 0 = a x 0          | 0 = -1 + 1
a x 0 = a x (-1 + 1)   | a x (b + c) = a x b + a x c
a x 0 = a x -1 + a x 1 | a x 1 = a
a x 0 = a x -1 + a     | -1 x a = -a and a x b = b x a 
a x 0 = -a + a         | -a + a = 0
a x 0 = 0

Now go forth and solve the problems of algebra!

Links for Tuesday

2010-09-08T00:33:00.003+03:00

A Blender-WebGL exporter got forked.

Rrrola's FX. Mandelboxes create a lot of cool geometry, riding that line between geometric and organic. I wonder if you could build game levels quickly by CSG cutting mandelboxes.

Some motion graphics videos on YouTube.

WebGL color spaces

2010-09-07T16:20:00.007+03:00

There's an interesting discussion happening on the public-webgl mailing list about the color spaces used by the textures and the canvas.

As far as I can tell, the problem is:

Simple math works correctly only in a linear color space.
Images and the rest of web content is 8 bits per channel sRGB.
Converting between two 8-bit color spaces is lossy.
You can't tell what an image's color space is from JavaScript.
You can't easily convert between color spaces in JavaScript.

The proposed solutions are:

Pass image data as-is to WebGL, treat WebGL canvas values as sRGB. (I.e. the do-nothing solution. If application writers want correct math, they should use linear textures and do a final linear-to-sRGB pass on the rendered image.)
Optionally convert textures to linear color space, treat WebGL canvas as linear. (Problematic as there is severe information loss at the conversions. Math works correctly though.)

Color spaces? But I thought space was black?

Well, actually, space has no color. It's transparent. It just looks black because there's very little light coming through it at most points.

Color spaces, on the other hand, are definitions of what the numerical pixel values in an image mean. A linear RGB color space means something like "given an output device with three component emitters at wavelengths R=580 nm, G=490 nm and B=440 nm, the number of photons emitted by an emitter is Max_Output * V, where V = component_value / (2^component_bits-1)." For sRGB, the component values are transformed by the sRGB curve (roughly pow(V, 2.2)) to have more control over the important parts of the dynamic range. In particular, sRGB sacrifices precision in the light end to have more control in the darks. Converting linear->sRGB loses detail in the light end, whereas sRGB->linear loses detail in the dark end.

To project an image onto your eyeball, the image data goes through the following wringer: Image data in image color space -> Viewer app color space -> Window system color space -> Display device color space -> Hardware blinkenlichts -> Eyeball.

Now, most software goes "huh?" at the mention of that and acts as if the image data is in the window system color space and uses the data directly. Further, the assumption when doing image editing operations is that the aforementioned color space is linear. (Why that might cause problems is that while 0.25 + 0.25 has the same luminance as 0.5 in a linear color space, to add the luminances in gamma 2.2 you have to do (0.25^2.2 + 0.25^2.2)^(1/2.2) = 0.34.)

Software that requires accurate color handling usually converts the 8-bit values to 16 or 32-bit linear color space to get around the conversion errors and to get easy correct math internally. The resulting high-depth image is then converted down to window system color space for viewing.

The real solution is to use 64 bits per channel and specify the channels to mean photon count / second / mm^2. With that, you can crank from 1 photon / second all the way up to 6500x noon sunlight. And linear math. Then use that as the interchange format between everything. Heyo, massive framebuffers~

In the meantime, maybe the browser could have an iconv-like system for converting images from one color space to another. WebGL 1.0 doesn't have float textures, so upconverting to linear float textures has to wait for the float texture extension.

With linear float textures, you could upconvert all images to float textures, then use a linear float FBO for rendering, and finally convert the FBO to sRGB to display on the canvas element.

Ramble on sorting

2010-09-06T23:09:00.006+03:00

Merge sort merges already sorted regions together. It also has an operation to turn a small region into a sorted region. The merge sort splits its input into small regions, turns each of those into a sorted region, and finally merges them together.

There are some sorted regions that have a trivial merge operation: if the largest value in region X is smaller or equal to the smallest value in region Y, the merge operation is the concatenation X ++ Y. We can pre-process the input of the merge sort so that the two merged regions have this property. The pre-processing pass uses a pivot value and moves all values smaller than the pivot to the start of the region, leaving the end of the region for values larger than the pivot. Now you don't have to do anything to merge the regions together once they're sorted. This is the idea behind the quicksort algorithm.

Because pivoting can cause the regions to have different sizes, there is a worst-case scenario where the pivot is the only element in the other region. Because of this scenario, quicksort has O(n^2) worst-case time complexity. If you pick the pivot to be the median of the region with an O(n) median finding algorithm and further split the equal-to-pivot values equally between the subregions, you can get an O(n lg n) quicksort (but the O(n) median finding algorithm has a high overhead, so they say the resulting algorithm is slower than quicksort in the average case).

You can do a O(n) pre-processing pass over the input data to find all already sorted regions in the input data. If the number of found regions is small, you're dealing with nearly sorted data. You can directly merge all sorted regions that are larger than the region size for the region sorting operation. The merging should be done so that the smallest regions are merged together first. Merging regions of unequal size causes the worst-time complexity to increase (merging a 2-region and a X-region takes X+2 time, consider N/4 2-regions merged into a N/2-region: it takes N/4 * (N/2 + (2+N) / 2) time, which is in O(N^2).)

The pre-processing pass can also reverse regions that are sorted in the opposite direction. And it is probably possible to topologically sort the sorted region descriptors by their min-max-values to find regions that can be merged by concatenation. Another merge trick would be to use binary search on one region to find the span that the other region overlaps, then you could get away with

memcpy(dst,smaller,smaller_sz); merge(dst+smaller_sz,overlap,region2,overlap_sz,region2_sz); memcpy(dst+smaller_sz+overlap_sz+region2_sz, larger, larger_sz)

.

If you keep track of how equally the quicksort pivot pass splits a region, you can avoid the quicksort worst-case by handing pathological splits over to merge sort. That is a bit nasty though, as the resulting regions no longer have the concatenative property and must be merged the slow way. Another alternative is to do something like introsort and switch over to heapsort once recursion depth exceeds lg n.

Once merge sort or quicksort reaches the maximum unsorted region size (say, 8 elements), implementations usually switch over to a low-overhead O(n^2) sorting algorithm like insertion sort.

The partitioning operation of merge sort is a no-op, but the partitioning operation in quicksort takes O(n) time. In the same fashion, the merge operation of quicksort is a no-op, but the merge operation of merge sort takes O(n) time.

You can do the merge sort merge in parallel by picking a pivot and splitting the sorted regions into two regions: one with values smaller or equal to the pivot and the other with values larger or equal to the pivot. The splitting operation consists of finding the index of the pivot in each region using a binary search. You then merge the small values with the small values and the large values with the large values. The merge between the merged smalls and larges is concatenation, so you can just allocate an array with the beginning reserved for the small values and the end reserved for the large values, and run the small and large merges in parallel. Recurse to increase parallelism. The nice thing about this approach is that you get segregated in-place parallelism, each of the sub-partitions can write to its own part of the root target array. The problem it suffers from is the quicksort partitioning problem: the pivot should be the median of the data set. Because the sorted regions are sorted, we can get the medians of each region easily enough, the worst case with one of those should be a 25:75-split when the regions are of equal size.

If your input data is 256 elements long and you split it into 8-element insertion-sorted regions, executing each in parallel (on a 32-core), what's the merge time? The first merge level you want to split into 2-way parallel merges (16 merges -> 32 merges). The second level is split 4-way (8 merges -> 32 merges), then 8-way and 16-way? The split bsearch overhead is 2 log₂ mergeSize per level (sub-splits execute in parallel).

To do the 16-way split (assuming we get a perfectly balanced split), we first do a 2 log₂ 128 = 14 op split, then 2 log₂ 64 = 12 op split, then 10 op and 8 op, for total split cost of 44 ops. The merge then takes 16 ops to run (8+8-merge), for total merge cost of 60 ops. If we stop the split one level short, it'd take 36 ops to split and 32 ops to merge, 68 ops in total. So I guess you do want to split it all the way down if you have negligible overhead.

The total runtime for the zero-overhead parallel merge tree would be 44 + 30 + 18 + 8 + 4*16 = 164 ops. Total sort time with 64-op 8-elem sorts would be 228 ops. Hey, sub-O(n). Quadruple the input size and the core count and it's 8+18+30+44+60+78+6*16 + 64 = 398 ops. Hey, sub-linear parallel scaling.

There is an another way to do parallel merging as well (cribbed from some lecture slides): for regions A and B, write A[i] to i + B.indexOfLarger(A[i]) and write B[i] to i + A.indexOfLarger(B[i]). N-parallelizable, serial O(N log N) time, parallel O(log N) time. However, it relies on scattered writes, so it doesn't seem very memory-access friendly (ideally each CPU core should be working on its own cache line, otherwise you get cache contention and traffic between the CPUs).

Fun reading: Heapsort, Quicksort, and Entropy.

Around the world on the ISS

2010-09-05T23:46:00.002+03:00

Two's complement integers

2010-09-04T05:50:00.007+03:00

Oh hey, I finally get how two's complement works after reading the first couple paragraphs of Why computers represent signed integers using two’s complement.

Simply put, adding a negative number to a positive number overflows if the result is positive, and doesn't overflow if the result is negative.

With 32-bit integers


5 - 1 = 5 + (-1) = 5 + (2³²-1) = 2³²+4 = 4
5 - 3 = 5 + (2³²-3) = 2³²+2 = 2
5 - 7 = 5 + (2³²-7) = 2³²-2 = -2

This is kinda interesting, we can put the separator between positive and negative numbers where we like:


void print_custom_int (unsigned int n, unsigned int separator)
{
  if (n < separator) {
    printf("%u\n", n);
  } else {
    printf("-%u\n", ~n+1);
  }
}

print_custom_int(0, 200);
// 0
print_custom_int(1, 200);
// 1
print_custom_int(-1, 200);
// -1
print_custom_int(200, 200);
// -4294967096
print_custom_int(-4294967097, 200);
// 199

Or maybe you need extra positive range and only 256 negative values:


print_custom_int(4000000000, 0xFFFFFF00);
// 4000000000
print_custom_int(2*2004001006, 0xFFFFFF00);
// 4008002012
print_custom_int(-32-80+15, 0xFFFFFF00);
// -97

// Watch out for overflows!
print_custom_int(-16*17, 0xFFFFFF00);
// 4294967024

And yeah, division doesn't work out of the box.


unsigned int div_custom_int (unsigned int a, unsigned int b, unsigned int separator)
{
  if (a < separator && b < separator)
    return a / b;
  else if (a < separator)
    return ~(a / (~b+1))+1;
  else if (b < separator)
    return ~((~a+1) / b)+1;
  else
    return (~a+1) / (~b+1);
}

unsigned int n20 = -20, n3 = -3, p4 = 4, p2 = 2;

print_custom_int(n20 / p4, 0xFFFFFF00);
// 1073741819
print_custom_int(p4 / n3, 0xFFFFFF00);
// 0
print_custom_int(p4 / p2, 0xFFFFFF00);
// 2
print_custom_int(n20 / n3, 0xFFFFFF00);
// 0

// vs.

print_custom_int(div_custom_int(n20, p4, 0xFFFFFF00), 0xFFFFFF00);
// -5
print_custom_int(div_custom_int(p4, n3, 0xFFFFFF00), 0xFFFFFF00);
// -1
print_custom_int(div_custom_int(p4, p2, 0xFFFFFF00), 0xFFFFFF00);
// 2
print_custom_int(div_custom_int(n20, n3, 0xFFFFFF00), 0xFFFFFF00);
// 6

The autumn is on its way

2010-09-02T01:51:00.002+03:00

Productivity tips, tricks and hacks for academics.

Branching factor generator

2010-08-31T23:31:00.003+03:00

A function to generate a stream of tree branching factors to approximate a non-integer branching factor.


branches' x i sum = 
  let v = (if sum/i > x then floor x else ceiling x) :: Int in 
  v : branches' x (i+1) (sum+fromIntegral v)

branches n = branches' n 0 0       

e_branches = branches (exp 1)
e_100000 = take 100000 $ e_branches
average l = sum l / fromIntegral (length l)
average $ map fromIntegral e_100000 == 2.71828

You could also write the average function as


average2 l = snd $ foldl (\(i, avg) e -> (i+1, (avg*(i-1) + e)/i)) (1, 0) l 

average2 $ map fromIntegral e_100000

even though that is less accurate. You can use a similar trick for the branches' function to use a running average instead of a running sum.


branches2' x i avg = 
  let v = (if avg > x then floor x else ceiling x) :: Int in 
  v : branches2' x (i+1) ((avg*(i-1)+fromIntegral v)/i)

branches2 n = branches2' n 1 0       
average $ map fromIntegral $ take 100000 $ branches2 (exp 1)

To generate a tree of wanted size out of the branching factor list, you'd have a tree generator that eats values off the list and generates tree nodes breadth-first. It is kind of a pain to write.

And some Tuesday links

2010-08-31T03:12:00.005+03:00

Man-Computer Symbiosis by J.C.R. Licklider. Continuing on the HCI vision golden oldies hitlist.

Twin hurricanes.

A Dust Mite.

Monday links

2010-08-30T20:27:00.002+03:00

Hawaiian lava flows look amazing via Highly Allochthonous.

Sunday links

2010-08-30T00:36:00.003+03:00

Even though it's technically Monday already.

Yellowstone: what lies beneath via Eruptions.

As We May Think by Vannevar Bush.

De architectura by Vitruvius.

The Periplus of the Erythraean Sea.

"My work in the houses of the Gods is abundant, and now I have begun an undertaking: Send much gold!"

Video of asteroid discovery, best viewed in high-def.

Hunza debris flow video.

Folding, part 3: parallel folds

2010-08-28T22:23:00.003+03:00

So, you can do a parallel fold if your folding function is associative.


-- Here i should be the neutral element of f, that is f i x = x.
-- E.g. 0 for + and 1 for *
parFold f i [] = i
parFold f i [x] = f i x
parFold f i xs =
  l `par` r `pseq` f l r
  where l = parFold f i left
        r = parFold f i right
        left, right = splitAt (length xs `div` 2) xs

-- See how associative functions work right with it
parFold (+) 0 [1..10] == foldl (+) 0 [1..10]
parFold (*) 1 [1..10] == foldl (*) 1 [1..10]

-- But others have problems
parFold (-) 0 [1..10] /= foldl (-) 0 [1..10]
parFold (/) 1 [1..10] /= foldl (/) 0 [1..10]

-- As - and / are defined as A + -B and A * 1/B,
-- we can work around these instances
parFold (+) 0 (map negate [1..10]) == foldl (-) 0 [1..10]
parFold (*) 1 (map recip [1..10]) == foldl (/) 1 [1..10]

But for the general case, we need something stronger. Roll in the reduceFold! It splits the foldable list into sublists, folds each sublist, and reduces the subresults into a single value.


reduceFold reduce f init [] = init
reduceFold reduce f init [x] = f init x
reduceFold reduce f init xs =
  l `par` r `pseq` reduce l r
  where l = reduceFold reduce f init left
        r = reduceFold reduce f init right
        (left, right) = splitAt (length xs `div` 2) xs

-- Associative operations work as usual
reduceFold (+) (+) 0 [1..10] == foldl (+) 0 [1..10]
reduceFold (*) (*) 1 [1..10] == foldl (*) 1 [1..10]

-- And we can now combine the subresults in a different fashion
reduceFold (+) (-) 0 [1..10] == foldl (-) 0 [1..10]
reduceFold (*) (/) 1 [1..10] == foldl (/) 1 [1..10]

Let's write some parallel list operations in terms of reduceFold.


parMap f lst = reduceFold (++) (\l i -> l ++ [f i]) [] lst
parFilter f lst = reduceFold (++) (\l i -> if f i then l++[i] else l) [] lst
parConcat lsts = reduceFold (++) (++) [] lsts
parConcatMap f lst = reduceFold (++) (\l i -> l++(f i)) [] lst

optLeft (Just x) _ = Just x
optLeft Nothing (Just x) = Just x
optLeft Nothing Nothing = Nothing

optFilter f i = if f i then Just i else Nothing

parFind f lst = reduceFold optLeft (\l i -> 
  optLeft l (optFilter f i)) Nothing lst

optRight x y = optLeft y x

parFindLast f lst = reduceFold optRight (\l i -> 
  optRight l (optFilter f i)) Nothing lst

-- Yeah, they pretty much work
divisibleBy n x = x `mod` n == 0

parMap (+ 10) [1..10] == [11..20]
parFilter (divisibleBy 3) [1..15] == [3,6,9,12,15]
parConcat [[1], [2,3,4], [5,6,7], [8,9], [10..20]] == [1..20]
parConcatMap (\x -> [x,x,x]) [1..3] == [1,1,1,2,2,2,3,3,3]
parFind (divisibleBy 127) [90,93..900] == Just 381
parFind (> 10) [1..10] == Nothing 
parFindLast (divisibleBy 127) [90,93..900] == Just 762

You might've noticed that both parFold and reduceFold are creating a parallel evaluation spark for each element in the lists. This is not a very good idea as the sparking overhead is often going to dominate the execution time and kill performance.

So we need some way to control the reduce tree branching factor and clump up the leaf-level folds into longer operations. I don't know if I can manage to come up with an automagic or even semi-automagic way to do this. I would like to use yesterday's theoretical math as guidance for the tree shape optimizer. Maybe we'll find out in the next installment.

A small hill-climb optimizer

2010-08-28T16:12:00.006+03:00


hillClimb f current delta error 0 = current
hillClimb f current delta error steps | delta < error = current
hillClimb f current delta error steps =
  if curVal > prevVal && curVal > nextVal
    then hillClimb f current (delta / 2) error (steps-1)
    else if prevVal > nextVal
           then hillClimb f (current-delta) delta error (steps-1)
           else hillClimb f (current+delta) delta error (steps-1)
  where curVal = f current
        prevVal = f (current-delta)
        nextVal = f (current+delta)

May have bugs in it, only tested on (\x -> -(x**2)).

Folding, part 2: parallel reduction math

2010-08-27T21:09:00.006+03:00

Parallel folding is a bit of an oxymoron. The fold operation sequentially builds up an accumulator value from the elements of a data structure. So if it's inherently sequential, how can it be parallelized?

If the reduce function given to fold is associative, the fold can be parallelized. If it's not associative, the fold is sequential.

Associativity means that you can turn (a + b) + c into a + (b + c), or, in function call terms, (f (f a b) c) == (f a (f b c)). When you have a long call chain like a + b + c + d + e + f + g + h, it means that you can split it into (a + b) + (c + d) + (e + f) + (g + h), evaluate the parenthesized bits in parallel, and then reduce the results to the final value. The reduction can also be parallelized, as it is the fold (a_b + c_d) + (e_f + g_h). This sort of parallel tree reduction is what map-reduce is all about btw: compute a bunch of values in parallel and reduce them to the result value.

Where it all goes by the wayside

Let's do the math to figure out the theoretical runtime for a reduction network. One parallel reduction step takes computeTime time. Moving the result to the next reduce node takes totalResultSize / bandwidth time. The variable totalResultSize is defined as resultSize * branchingFactor, as the branching factor is the amount of nodes sending their reduce results to the reduce node. We can also split it into a more useful form branchingFactor * (resultSize / bandwidth) and write transferTime = resultSize / bandwidth Finally, the amount of parallel reduction steps required to reduce a tree of size N is log_branchingFactor N (and 1 < branchingFactor <= N).

Putting it all together, we get something like


treeReduceTime n branchingFactor computeTime transferTime =
  stepCount * stepTime
  where stepTime = computeTime + (branchingFactor * transferTime)
        stepCount = if branchingFactor == 1.0
                      then n
                      else log n / log branchingFactor

From the above we can come up with a couple of generalizations:

If computeTime is O(1) with regard to branchingFactor and you have infinite bandwidth or zero result size (both mean that transferTime is zero), you want a branching factor of N. That's easy enough to see:


With infinite bandwidth, the treeReduceTime is
  (computeTime + (branchingFactor * 0)) + stepCount =
  computeTime * stepCount

To minimize the above, we need to minimize stepCount.
The smallest stepCount we can have is 1.
The stepCount is 1 when branchingFactor is N.

If you have infinitely fast compute nodes (computeTime = 0), you want a branching factor of e. That's a bit tougher to prove[1] (or maybe I'm just lacking some essential bit of knowledge that would make the proof easier.)

Note that it's highly likely that computeTime depends on branchingFactor. If the reduce algorithm is O(n), computeTime = branchingFactor * elementComputeTime, and in that case the optimum branching factor is e (as elementComputeTime and transferTime have the same factor, they can be rolled together and the infinitely fast compute nodes proof applies). For O(n^2) algorithms, the optimal branching factor seems to be sqrt(e).

For O(1) computeTime, finite computational power and non-zero transfer time, the optimal branching factor is somewhere between e and N. I used a simple hill-climb optimizer to find optimal branching factors for different computation/transfer ratios, and they do seem to be independent of N (well, apart from the maximum branching factor). For 1:1 compute time : transfer time, the branching factor is around 3.6. For 1:10 (low bandwidth), the branching factor is a bit over 2.8. For 10:1 (slow compute), the branching factor is about 8.6. Here's some approximate ratios for integer branching factors: 3 -> 1:3.5, 4 -> 1.5:1, 5 -> 3:1, 6 -> 5:1, and for 1000:1 the factor is about 226.

The above approximations can be used for branching factor -dependent computation times too. You do that by replacing the computation time with per-reduce-step overhead and combining computeTime in transferTime. The new ratio is overhead : (computeTime+transferTime). For intuition on this: If the overhead dominates, you want more branching as that will decrease the contribution of the overhead. If you have very little overhead, you want a branching factor close to e.

For an example, consider an email tournament that aims to pick the best email out of ten thousand. You receive some fixed amount of emails and your task is to pick the best one, then forward it to your assigned receiver (who does the same, etc., and the last reducer sends the mail to the award ceremony). Let's estimate that the transfer time is 5 seconds, the computing time is a minute per email, and the fixed overhead is 8 hours. This gives us an overhead:processing-ratio of 443. If each receiver processes 10 emails per session, running the tournament would take 33 hours. But if each receiver does 118 emails a go, the tournament would be done in 19.5 hours. And if they do 200 emails, it'd take ... 20.2 hours. (Yeah, 118 is the best branching factor here.)

Suppose you do the same tournament but with people glued to their mail apps. Now the overhead might be just one minute. This time, doing the tournament with 118 mails per session would take a snappy 4.1 hours. But as the role of overhead has fallen, we should reduce the branching factor as well. With 4 mails per session, the tournament would take only 35 minutes. (With 2 mails per session, 42 minutes. 4 is optimal.)


With infinitely fast compute nodes, the treeReduceTime is
  (branchingFactor * transferTime) * stepCount

For a branching factor of 1, that equals
  transferTime * n

For other branching factors, the equation is
  b * t * log_b(n), where b = branchingFactor and t = transferTime

From which we see that t is a fixed factor,
so the minimum of the above equation is the minimum of the equation
  b * log_b(n)

Which can be written as
  log_b n = log n / log b,  b > 1
  b * log_b(n) =
  b * (log n / log b) =
  (b / log b) * log n

From which we see that log n is a fixed factor,
so the minimum of b * t * log_b(n) is the minimum of
  b / log b

To find the minimum, find the derivative
  D(b/log b) = (1*log b - b*1/b) / (log b)^2
             = (log b - 1) / (log b)^2

And its inflection point at b = e,
  (log e - 1) / (log e)^2 = (1 - 1) / 1^2 = 0 / 1 = 0,
we check the values of b / log b on both sides of the
inflection point to be larger than the value at the inflection point,
  e / log e = e =~ 2.718
  2 / log 2 =~ 2.885
  3 / log 3 =~ 2.731,
thus confirming it as a local minimum.

Checking the discontinuity of b / log b at b = 1,
  b / 0 = infinity,
and discovering it to be larger than e / log e,
we find the global minimum of b / log b to lie at b = e.

Plugging that back into the original equation:
  minimum of branchingFactor * transferTime * stepCount is
    transferTime * min(n, e*ln(n))
  and because n >= e*ln(n), the minimum is
    transferTime * e * ln(n).

Folding, part 1

2010-08-26T16:03:00.003+03:00

The fold function is a data structure accumulator: it goes through every element in the data structure and builds up an accumulator value from them. An example for folding lists:


foldl f acc [] = acc
foldl f acc (x:xs) = foldl f (f acc x) xs

With a fold you can do maps and filters of different kinds. Here's a few examples:


-- id x = x
reverseMap f lst = foldl (\l i -> (f i):l) [] lst
reverse = reverseMap id
map f lst = reverse (reverseMap f lst)

filter f lst = 
  reverse (foldl (\l i -> 
    if f i 
      then i:l 
      else l
  ) [] lst)

concat a b = foldl (\l i -> i:l) b (reverse a)

concatMap f lst = foldl (\l i -> concat l (f i)) [] lst

triplicate lst = concatMap (\i -> [i,i,i]) lst

drop n lst = 
  reverse $ snd $ foldl (\(c,l) i ->
    if c >= n 
      then (c, i:l) 
      else (c+1, l)) (0,[]) lst

You can also do find operations with fold, but they're not as efficient as you'd like, because the fold goes through all the elements in the list. To make find efficient, we need a fold with break.


findInefficient f lst = 
  foldl (\e i -> 
    case e of 
      Nothing -> (if f i then Just i else e)
      Just x -> e) Nothing lst

data Break a = Break a | Continue a

foldlWithBreak f acc [] = acc
foldlWithBreak f acc (x:xs) = 
  case f acc x of
    Break a -> a
    Continue a -> foldlWithBreak f a xs

find f lst = 
  foldlWithBreak (\e i -> 
    if f i 
      then Break (Just i) 
      else Continue Nothing) Nothing lst

take n lst = 
  reverse $ snd $ foldlWithBreak (\(c, l) i -> 
    if c >= n 
      then Break (c,l)
      else Continue (c+1, i:l)) (0, []) lst

any f lst = 
  -- maybe False (const True) (find f lst)
  case find f lst of
    Just x -> True
    Nothing -> False

all f lst = 
  -- not $ any (not.f) lst
  not (any (\x -> not (f x)) lst)

And that's it for today. Tomorrow, parallel folding.

FormData file uploads

2010-08-26T01:14:00.005+03:00

XMLHttpRequest Level 2 FormData is super. Loving it.


var fd = new FormData();
fd.append("key", "value");
var files = dropEvent.dataTransfer.files;
for (var i=0; i<files.length; i++)
  fd.append("file", files[i]);

var xhr = new XMLHttpRequest();
xhr.open("POST", "/upload");
xhr.send(fd);

sse.h update, HTML5 fiddling, more links

2010-08-25T02:31:00.002+03:00

I updated sse.h with license info (MIT) and added software implementations of recip() and rsqrt() to double2 and double4. This Vector Classes library might be more useful for Real Serious Use, though.

Currently fiddling with HTML5 media tags (cursed be audio and video codecs. If it were possible to implement media codecs in JS, that would be a solution.) Also doing drag-and-drop file uploads. The event API is ... interesting. To receive a drop event, you need to preventDefault other drag-related events.

Failures are strengths.

Have some links to go:

Eruptions blog readership.
Activity at Sakurajima Volcano Intensifies.
Region-based memory management, ML Kit, A Simplified Account of Region Inference, Combining Region Inference and Garbage Collection, A Safe Runtime System in Cyclone.

A new blog post!

2010-08-23T23:43:00.004+03:00

No, wait, sorry, just more link copypasta. My bad.

Join, or Die. Benjamin Franklin, truly a master of Perl. Via Hark, a vagrant.

Ice Island Calves off Petermann Glacier, via Bad Astronomy.

Yet more random links

2010-08-22T23:25:00.004+03:00

Electric blue planktom blooms off Ireland.

Employer Expectations, Peer Effects and Productivity: Evidence from a Series of Field Experiments [PDF]. A study using Amazon Mechanical Turk, one finding was that peers reward effort rather than output quality. Quote from the implications section: "For example, it may be difficult to get workers to substitute easy, correct procedures for difficult, inefficient procedures. Ironically, the difficulty itself might make an outdated procedure harder to replace, as workers who adopt the easier method might be perceived to be shirking. The finding that exposure to low-output work lowers output, combined with the finding that low-productivity reduces willingness to punish, suggests the possibility of an organizational vicious cycle: after observing idiosyncratically bad work, workers may lower their own output and punish less in response, in turn reducing other workers’ incentives to be highly productive."

Forecasting bird migration for the Finnish Air Force. "Why? Bird is a threat."

Random links for another day

2010-08-22T01:47:00.004+03:00

PacketShader: a GPU-Accelerated Software Router - eeh?

Enceladus plumes with anigifs. Cassini just keeps taking these awesome shots of the Saturn system.

Kirkenes harbour (and googling with the name of the ship got me this. Take that for good measure as well.)

Fantoft stave church, Pre-Colombian Mexican architecture

... the deep-sea food web is fueled by a rain of dead plants and animals from surface waters. Sediment for the sediment god, let the red clay sink back into the mantle.

Cryosphere Today.

That reminds me: We inhabit a few meter thick layer at the solid-gas interface of a 6500 km radius sphere of rock and iron. Go up or down a few meters and you'll die. Go to the solid-liquid or gas-liquid interface and you'll die. Go to the solid, go to the liquid or go to the gas and you'll die.

Bismuth crystal. I quite like the Mineralogical Record print mag. Lots of good pictures.

Random links of the day

2010-08-21T01:02:00.007+03:00

Use of Ada in a Student CubeSat Project (via Ada-Europe 2008 via SPARK Conference Papers). First they built some arctic sea ice buoys and then used the same HW/SW-platform to do a CubeSat, in a high-integrity dialect of Ada. (The conference proceedings have some other wacky stuff too, like a presentation about porting a naval command system used in Royal Navy destroyers to Ada 2005. And yes, the slides are basically a long list of stuff that's changed between language versions, so it's not that exciting.)

Neptune about to complete first orbit since discovery. Discovered in 1846, orbital period 165 years, orbit complete next summer. Celebration? (Via Uncertain Principles.)

The dangers of plein air. That reminds me, I've fallen off the drawing routine. Perhaps I should unfall. Also this.

"According to the meteorologists, it's never summer in Gamvik". The coast of the Arctic Sea is a jolly good place. Full of fun and excitement and days that last for months.

Random programming ideas

2010-08-19T13:13:00.006+03:00

Programming is about gathering knowledge about a piece of code. What does it do, what are its limitations, how does it perform, where can it be used, what is it doing right now, how do I extend it, how much memory does it use, how do I know that it works right, etc.
Testing, logging, documentation and profiling as an integral part of a language. Parse each function and module into a {body, tests, docs, logstream, profilestream}-struct that can be fiddled with at runtime.
Automatically test a function based on its type. Call it with different args, log the results and summarize detected patterns to give programmer a quick description of what the function does and what its time-space-complexity looks like. Add a way to flag found properties as want/do not want and check that subsequent versions of the function fulfill those (basically autogenerating QuickCheck rules).
Typed UNIX pipes. Use HTTP headers to tell the recipient what's coming through.
Auto-generate type->type-pipelines by putting conversion functions in a graph, run graph search on it with edge weight given by amount of information lost in each conversion (and secondarily performance).
Semantic type information to encode the effects of a function. Two functions that purport to do the same thing may well have semantic differences: suppose two image scaling functions where one assumes linear color space and the other does gamma correction. Both do ((RGB w h), new_width, new_height) -> (RGB new_width new_height), but the results differ.
Concurrency based on HW-level message passing. Processors read and write in cache lines, cache lines can't be shared between two processors => use cache lines as messages.
Array programming-like parallel array computation lib. Element-wise operations on N-wide vectors are easy to parallelize and vectorize: divide N by the amount of processors and the SIMD width. Needs to compile array-level pipes ((f . g . h) A) into element-level pipes (map (f' . g' . h') A) to maximize arithmetic per memory access. Maybe Data Parallel Haskell does this stuff already.
Uniform batteries-included interface for different data structures. Yeah, doable with typeclasses. Yeah, something like OCaml Batteries Included Enums. You can reduce data structures to either empty, fold and unfold (for streaming data structures) or alloc, length, get and set (for random access structures).
Data structure -generic regular expressions (i.e. not tied to just strings), are they useful? See proof of concept. What you need is a fold over the data structure (and I guess that most data structures are foldable, how else would you free the memory their elements use?)
Instruction set with bounded & bordered LOADs and STOREs: LOAD reg addr minAddr maxAddr [borderValue]. Throw hardware exception if a bounded LOAD goes out of bounds, load borderValue if a bordered LOAD goes out of bounds. Bordered STORE writes to the address in borderValue instead. Faster than two compares and JMP? More compiler-friendly? Eliminates buffer overflows? (Random literature search: Security improvement in embedded systems via an efficient hardware bound checking architecture, Checking Array Bound Violation Using Segmentation Hardware, A Comprehensive Approach to Array Bounds Check Elimination for Java, x86 segmented access mode, some mainframe architectures are said to have something like that, mysteries of life.)

China at #2

2010-08-16T06:08:00.007+03:00

Hey, PRC thundered past Japan in nominal GDP this year. At this pace, China'll be passing USA and EU in a decade and maybe reach the same GDP per capita by 2040. At which point its annual growth rate will stabilize at around 4% and all the funds expecting 10% will cause a major economic crash, absorbed in part by an ascendant India? I guess the question is how China will transform from a fast-growing implementation economy into a slow-growing research one (IN THE FUTURE!!!1!one).

A slow cache algorithm and a faster one

2010-08-10T01:33:00.003+03:00

It took around 30 ms to draw some text titles. The weird thing was that not drawing the titles didn't really affect the draw time. So clearly the problem was not in the drawing speed. Investigating, my text bitmap LRU cache was hitting its max size, destroying all possible benefits of caching. But making the cache large enough didn't fix the performance problem. A mystery!

On closer look, my caching algorithm was retarded:


function makeCache(maxSize) {
  return { array: [], hash: {}, maxSize: maxSize };
}

function cacheItem(cache, item) {
  // shift old items off the array
  while (cache.array.length >= cache.maxSize) {
    var it = cache.array.shift();
    delete cache.hash[it.key];
  }
  cache.array.push(item);
  cache.hash[item.key] = item;
}

// move item with key to the end of cache array
function refresh(cache, key) {
  var it = cache.hash[key];
  var idx = cache.array.indexOf(it);
  cache.array.splice(idx,1);
  cache.array.push(it);
}

function get(cache, key) {
  refresh(cache, key);
  return cache.hash[key];
}

Why that is bad: O(n) get, O(n) cacheItem after hitting maxSize. If you're going through the whole cache every frame, that's O(n^2). If you're also at maxSize, double that! Aargh!

So I rewrote it to be less insane and now the drawing is nice and fast:


function makeCache(maxSize) {
  return { array: [], hash: {}, maxSize: maxSize, lastUsedIndex: 0 };
}

function cacheItem(cache, item) {
  // chop off the oldest items after reaching overflow limit
  if (cache.array.length >= cache.maxSize * 2) {
    cache.array.sort(function(a,b){ return a.lastUsed - b.lastUsed; });
    var deleted = cache.array.splice(cache.maxSize);
    for (var i=0; i<deleted.length; i++) {
      var it = deleted[i];
      delete cache.hash[it.key];
    }
  }
  cache.array.push(item);
  cache.hash[item.key] = item;
}

// update the cached item's lastUsed
function refresh(cache, key) {
  cache.hash[key].lastUsed = cache.lastUsedIndex++;
}

function get(cache, key) {
  refresh(cache, key);
  return cache.hash[key];
}

The gets are now O(1) and cacheItem with maxSize cache runs in amortized O((2n log 2n) / n) time. Or something. You could further optimize it by using a split instead of sort. Basically run an O(n) selection algorithm to find the maxSizeth element of the cache array, then use that as a pivot and do a single quicksort partition pass. That'd give you amortized O(2n / n), which is close enough to O(1) :P

Berlin!

2010-07-27T00:15:00.006+03:00

Hanging out in Berlin this week, with a shiny new sketchbook to boot.

Went to the Pergamonmuseum on Sunday to be all Ozymandias. The things they have there are huge (indoor display of 15m high pillars, etc.) and pretty old (2-3 thousand years) and beautiful. And messed up in that ancient style of messed-upness. The objects in the Islamic museum side were nicee, superior metalcraft and calligraphy (also a millennium or two younger than the rest of the stuff but who's counting). The Assyrian reliefs had writing going across them, like a textured strip of cuneiform. They had a special display on the coloring of the Greek statues, with brightly colored reproductions of much RGB might.

Silly C hacks

2010-07-20T22:27:00.007+03:00

Macros! To do safe mallocs and fopens! Brittle like the glass you tread upon!


/* A small cat that prints out the files given to it */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <malloc.h>
#include <sys/stat.h>

#define BLOCK(pre, post, ok, error, block) \
  { \
    pre; \
    if (ok) { \
      block; \
      post; \
    } else { \
      error; \
    } \
  }

#define with_malloc(t, var, sz, block) \
  BLOCK(t var = (t)malloc(sz), free(var), var, \
        fprintf(stderr, "Malloc failed, out of memory\n"); exit(1), block)

#define with_file(var, fn, mode, error, block) \
  BLOCK(FILE* var = fopen(fn, mode), fclose(var), var, error, block)

#define with_file_contents(var, len, fn, error, block) \
  with_file(var##__f, fn, "r", error, \
    off_t len = file_size(fn); \
    with_malloc(char*, var, len+1, \
      len = fread(var, 1, len, var##__f); \
      var[len] = 0; \
      block))

off_t file_size(char *path)
{
  struct stat st;
  stat(path, &st);
  return st.st_size;
}

int main(int argc, char *argv[]) {
  int i;
  for (i=1; i<argc; i++) {
    with_file_contents(s, slen, argv[i], fprintf(stderr, "Couldn't open %s.\n", argv[i]),
      fwrite(s, 1, slen, stdout);
    );
  }
}

It would be kinda fun to port the whole functional array-munging part of prelude.ml to C. Fun in a very painful fashion.

Small C++ SSE3 vector lib

2010-07-20T18:21:00.005+03:00

If you want to do some simple and somewhat fast vector math in C++, sse.h might be fun.

It has float4, double2 and double4 structs, constructors of different kinds, operator overloading and methods for horizontal add, dot product and swizzling.

Code examples:

4x4 matrix multiply with floats and doubles:


// dst = a x b, ~44 cycles with -O3 -funroll-loops
void mmul4x4 (const float *a, const float *b, float *dst)
{
  for (int i=0; i<16; i+=4) {
    float4 row = float4(a) * float4(b[i]); // float4(a) is {a[0], a[1], a[2], a[3]}
                                           // float4(b[i]) is {b[i], b[i], b[i], b[i]}
    for (int j=1; j<4; j++)
      row += float4(a+j*4) * float4(b[i+j]);
    *(float4*)(&dst[i]) = row;
  }
}

void mmul4x4d (const double *a, const double *b, double *dst)
{
  for (int i=0; i<16; i+=4) {
    double4 row = double4(a) * double4(b[i]);
    for (int j=1; j<4; j++)
      row += double4(a+j*4) * double4(b[i+j]);
    *(double4*)(&dst[i]) = row;
  }
}

Dot product of two arrays of vectors:


double dotArrays (const double2 *a, const double2 *b, int len)
{
  double2 sum;
  for (int i=0; i<len; i++) {
    sum += a[i] * b[i];
  }
  return sum.sum();
}

AC the planet

2010-07-12T16:46:00.004+03:00

+7C equals ow. (Check out the map in that: yellow areas would be deadly.)

Guess humans wouldn't have survived too well in the Cretaceous period (when the avg. temp was 15-20C higher than today. Tropical sea surface temperatures at around +37C. Here a fun lecture on that.) If you'd time-travel back to ride on brontosaurs or whatever, you'd probably end up with a heat stroke.

Future display tech

2010-07-09T23:13:00.008+03:00

Graphics, where is the limit? The display resolution of the day is 1920x1080, or Full HD, as it is called, with high-end displays operating at 2560x1600. Modern graphics cards are quite capable of running games at those resolutions. Adding stereo-3D on top of that doubles the amount of pixels, but that isn't much of a problem either. Pushing up the resolution from 100 ppi to 300 ppi would need an extra 9x the resolution, which would probably start pushing today's limits.

To do no-glasses viewing-angle independent 3D requires a pixel for every direction (or something like a pixel, e.g. plop a lens in front of a flat pixel array). That would easily multiply the resolution by 16 or 256 or such. Not to mention requiring new display tech and the cameras to produce content for them.

You probably want to push the dynamic range of the display as well, in order to reproduce bright noonday sun and a pitch-black night. According to this page the contrast from the sun to the limits of human night vision is around 300 trillion to one. To do that range with fixed step would take 48 bits. The display would have to be able to emit from a couple photons at a time all the way up to a billion candela per square meter, which is a good million times brighter than today's bright LCD displays.

Using 16 bits per RGB channel for encoding the hue and saturation, and then 48 bits divided among them to do the luminosity would result in 96 bits per pixel. With 32-bit alpha channel, 128 bits per pixel, or four times today's display bit depth.

The display would be nice if it covered your field of vision from a comfortable viewing distance. A wall-sized display might be the ticket to that. Let's say 5 meters by 3 meters.

Summing up all the above, a future wall-sized 256x integral hologram display with 300 ppi resolution and natural dynamic range would have 550 gigapixels at 128-bit color, and would need a 8.8 TB frame buffer from the graphics processor. A more modest 30" display at those specs would have only 9.4 gigapixels and a 150 GB frame buffer. If we assume that 1 GFLOPS is enough to run a 32-bit two megapixel display, a 128-bit 550 gigapixel display would require around 1,100,000 GFLOPS of processing power. Today's graphics cards are at around 2 GFLOPS, so if they keep doubling in performance every year it'd take 20 years for them to get to the required 1.1 PFLOPS.

What would such a display look like in person? I suppose the experience would be like looking out of a window. A window to a place that doesn't exist. Where the searing hot twin suns rise above the desert, windswept dialog boxes of long-forgotten programs casting long shadows into your living room. Or the view out into a starlit forest: a ghastly green shape here, another there, the darkness registering as a noisy shimmer as the rods in your eyes struggle to resolve the scant errant photons reflected from your surroundings.

I don't really know if 256x holography would be enough, though. Maybe you want to boost that to 4096x or more. And make the display wrap around you for that VR cave experience. And add one extra color channel for tetrachromats. Or heck, just use photon energy histograms as your colors instead of the limited RGB. The people with those new-fangled omni-EM eyes will thank you for that. And a monitor just isn't a monitor unless you can do astronomical spectrography off it, take your X-ray with it, microwave an egg, make your radio play tunes and do a realistic documentary of the Chernobyl disaster, complete with the viewer receiving a half dozen gray of radiation.

Optimization madness

2010-06-27T06:51:00.093+03:00

This paper, Believe it or Not! Multicore CPUs can Match GPUs for FLOP-intensive Applications!, from IBM Research is kinda fun. They did a benchmark with an Nvidia GTX 285 pitted against a dual quad-core Nehalem Xeon and a Power7. The benchmark is doing O(n²) dot products from two 4 MB (not MiB, mind) source images to a 250 kB destination image. Maybe you can see where this is going?

TL;DR

The aggregate cache bandwidth of two high-end $1386 quad-core Xeons almost reaches the main memory bandwidth of a $300 graphics card. The cache bandwidth of an eight-core Power7 blade (cheapest system I could find started at $9788.00) actually exceeds the main memory bandwidth of the graphics card by 40%. A cache-conscious version of the GPU implementation might make it 5x faster (as according to this blog post, the GTX 285 has 1.5 TB/s cache bandwidth). Currently trying to write one. Testing on my 2GHz Core 2 Duo and a Radeon HD 4850 that has 40% of GTX 285's bandwidth.

My results as of now

Hardware tested

ATI Radeon HD 4850, 67.9 GBps memory bandwidth, 384 GBps L2, 480 GBps L1, 1000 mul-add SP GFLOPS, hobbled by lack of local mem and the lack of image support in AMD's OpenCL drivers, resulting in most reads coming from main memory
Core 2 Duo E6400, 6.4 GBps memory bandwidth, 18 GBps L2, 90.5 GBps L1, 37.5 mul-add SP GFLOPS, hobbled by lack of cache bandwidth and computing performance (but hey, maybe there's some way to unlock the last 16 GFLOPS)

Table of current results (8 FLOPS per 32 bytes memory read)


  The results from the paper are on 500x500 input images. The paper reported
  runtimes in seconds, I calculated the bandwidth numbers by dividing 281 GB
  by the reported runtimes.

  The OpenCL results are measured without the context and kernel creation
  and release times, as those are likely only relevant for a single-run
  workload (and a 0.7 sec one-time overhead to do a 1 sec operation would
  be a bit skewy). The OpenCL results do include the data write and read times.
  GLSL times are similarly measured without global setup overhead.

  359 GBps (90 GFLOPS) 2x2-blocked GLSL on Radeon HD 4850 (1024x1024 input)
  337 GBps (84 GFLOPS) 2x2-blocked GLSL on Radeon HD 4850 (500x500 input)
  278 GBps (70 GFLOPS) optimized OpenCL GPU on Radeon HD 4850 (640x640 input)
 [276 GBps (69 GFLOPS) the paper's optimized impl on 8-core 3 GHz Power7]
  259 GBps (65 GFLOPS) optimized OpenCL GPU on Radeon HD 4850 (512x512 input)
  243 GBps (61 GFLOPS) optimized OpenCL GPU on Radeon HD 4850 (500x500 input)
  240 GBps (60 GFLOPS) 2x2-blocked GLSL on Radeon HD 4850 (250x250 input)
 [230 GBps (58 GFLOPS) the paper's GTX 285 implementation using local mem]
  204 GBps (51 GFLOPS) naive GLSL implementation on Radeon HD 4850
 [159 GBps (40 GFLOPS) the paper's GTX 285 implementation]
 [150 GBps (38 GFLOPS) the paper's optimized impl on 2x4-core 2.93 GHz Nehalem]
   86 GBps (22 GFLOPS) naive OpenCL GPU implementation (1D workgroup)
   84 GBps (21 GFLOPS) optimized OpenCL CPU on Core 2 Duo E6400
   45 GBps (11 GFLOPS) manually optimized OpenMP SSE implementation
   25 GBps (6 GFLOPS) naive OpenCL GPU implementation (2D workgroup)
    5 GBps (1 GFLOPS) naive OpenMP SSE implementation
    5 GBps (1 GFLOPS) naive scalar C implementation
    1.1 GBps (0.3 GFLOPS) naive OCaml implementation
    0.03 GBps (0.008 GFLOPS) naive pure-Python implementation

Notes

The GPU data transfer overhead scales linearly with problem size, but the computational complexity of the algorithm grows quadratically. So on small problem sizes the effect of the overhead is more visible, as you can see above. The CPU does not suffer from the data transfer and setup overheads, so there is a point where the CPU is faster for problems smaller than that. With the current results for the GLSL implementation and the server CPUs, that point is at somewhere around 200x200 input image size. For my Core 2, it's around 150x150. For single image runs (with all the OpenGL overhead factored in the runtime), the GLSL-Power7-cutoff runs at around 600x600.

CPU performance went up from a 5.5 GBps naive parallelized SSE implementation to a 84 GBps OpenCL optimized version (93% of L1 bandwidth, 56% peak GFLOPS). The fastest manual SSE code I managed topped out at 43 GBps. GPU performance went from a 86 GBps naive implementation to a 337 GBps optimized version (70% of L1 bandwidth, 8% peak GFLOPS). Note that the naive CPU implementation is written in C, parallelized and vectorized. If you're writing it in pure Python instead, the performance is around 0.03 GBps, or 180x slower than the C baseline. For another baseline, a naive OCaml implementation manages 1.1 GBps. And the OCaml and Python implementations don't get faster from reordering the memory accesses because they're bottlenecked by computation, not memory access.

But in vectorized C, this algorithm is all about optimizing your memory accesses. The CPU got a 6x boost from changing the loop nesting (4x from the order optimization, rest from it enabling each core use its L1 caches and thus making parallelization more effective) and another 2x boost from reusing already loaded values (again, reducing average read latency). The GPU seems to be mostly bound by the memory access pattern (and the lack of local memory and global mem caching on the Radeon 4xxx series). On the GPU I got a 2.2x performance increase from updating 8 values at a time and reading the current block of used values to a local array. Heck, a more optimal OpenCL access pattern with _more_ reads and computation (but, sadly, _wrong_ reads) nearly doubles the GPU performance. A mostly lock-step-iterating version that does more reads got me an extra 20 GBps on OpenCL (kinda freaky, that.) GLSL does use the caches, so it's a good bit faster.

The AMD OpenCL implementation in Stream SDK 2.1 doesn't support images on HD 4850, and global memory accesses aren't cached (if I understood correctly. The generated assembly has texture fetches with the UNCACHED keyword, so probably!), which might be the cause for the low bandwidth. Guess I should buy an HD 5xxx / Nvidia and test the performance with images (I hear they have real local memory as well :P). Testing with a faster CPU would be nice as well.. Not that I have money for any of those so pfft.

Some interesting things from my POV: Slow CPU code is slow. It's possibly to reach at least 93% L1 bandwidth on the CPU, though even that doesn't max out the ALU (56% of the peak GFLOPS). And it's still kinda slow even when you do reach it. OpenCL is pretty nice for writing high-performance CPU code, it takes the pain out of SSE and parallelization. The setup is a pain though, and there's a 0.1 s binary program loading overhead. GPUs are kinda fast, even when hobbled with incomplete driver support. I'm stuck at 8% the peak GFLOPS because I can't read from the memory fast enough (and the GLSL implementation doesn't reuse loaded values to compute several output values in parallel, so it's doing cache reads instead of register reads. Which might well help.)

The algorithm used in the benchmark is massively amenable to optimization as it has an embarrassingly parallel map pass, a completely associative reduce pass, completely vectorizable, each output value is basically a single N-wide mul with a horizontal add, has a large amount of shared reads between adjacent output values, benefits a lot from cache, benefits a lot from prefetch, etc. Any numbers you see about a benchmark like that should be taken with a huge grain of salt, it's just too easy to get 2x performance differences with small optimizations. The difference between my scalar implementation and parallelized SSE implementation was 15%. The difference between the parallelized SSE implementation and the optimized implementation was 1700%. The GPU performance delta was from a pessimized version that ran at 20 GBps to a regular version that ran at 85 GBps to an optimized version that ran at 240 GBps to a GLSL version that ran at 337 GBps. And the OpenCL GPU number is very sensitive to data access pattern, problem size and the size of the OpenCL workgroup: I'm getting 240 GBps with 64 workers on 512x512 images, 190 GBps with 32 workers, 180 GBps with 128 workers, 30 GBps on 500x500 images (that rises to 150 GBps by changing stride size), 150 GBps on 256x256 images, etc. The GLSL version is much stabler with regard to input size, though it too suffers from lower performance at smaller input sizes.

The CPU is less prone to wild swings in throughput from small changes, which is nice. The only way to get decent performance on this algorithm seems to be memory access optimization, SSE, parallelization and collecting several values at a time. So you pretty much need to write it in C / C++ with OpenMP and SSE intrinsics, or OpenCL. And if you're writing it in OpenCL, why not get the GPU in the mix as well...

GPUs would be so much better if they didn't hard-lock when encountering an infinite loop. Over the past few days I've done Alt-SysRq-REISUB more times than I can remember. The OpenCL GPU implementation on the HD 48xx is ... finicky. Small changes in code can lead to either the GPU hard-locking, a 10x performance drop, or just plain incorrect behaviour. I'd like to see either GPUs with CPU levels of usability or CPUs with GPU levels of performance. Or both. I'm not picky.

Based on my initial experience, I'd say that if your algorithm is amenable to GPUs (i.e. it does lots of single-precision FP math or needs high memory bandwidth), it's not too difficult to have a naive GPU implementation that performs better than a naive CPU implementation, especially if your CPU implementation is in a slow language. In my case, optimizing the CPU and GPU implementations happened in parallel and the optimizations were quite often portable between the different implementations. The CPU was easier to optimize for and had a much better development experience (i.e. no lockups), but my cheapo CPU didn't have the performance to match the GPU at any point. The CPU-GPU performance delta went from 15x between naive implementations (40x with GLSL) to 4x between the optimized implementations.

In terms of price-performance, a $100 GPU and a $100 quad-core CPU (hmm, should you include the memory price in the CPU prices?) should perform about the same on this algorithm when using OpenCL. With GLSL, the GPU would be around 50% faster. A 2.9 GHz Athlon II X4 should get something like 220 GBps (247 GBps L1 bandwidth) and a 3 GHz quad-core Core i7 should get around 175 GBps (192 GBps L1 bandwidth). At the $200 price point, the gap between the GPU and CPU will likely increase: you get around 2x the GPU L1 bandwidth, but only a 50% extra CPU L1 bandwidth (with 6-core Phenom. Opting for a higher-clocked quad-core CPU would likely show very little performance benefit, you'd only get +20% out of going from 2.9GHz to 3.5GHz, for example.) But lacking the hardware to test, this is mere speculation.

Development log

Update #0: Source code.

Update #1: GPU at 86 GBps. Managed to drop the CPU runtime from 49 s (5.5 GBps) to 30 s (9 GBps) by making the inner loop update the values in blocks. The intuition there being that you're going to have better cache locality if you're updating a bunch of values that use data from the same area. It also returns bad values for some input sizes, so it probably doesn't work right.

Update #2: The OpenCL GPU kernel is kinda hard to optimize, most of my naive tries ended up hurting performance. But I did manage to get the performance up by 40% using a similar block-based approach as on the CPU (by using a local cache array instead of relying implicitly on L1.) Currently at 118 GBps with roughly the same optimization effort as with the CPU implementation. I don't really know much anything about optimizing OpenCL though, so I think the number could go considerably higher. On the CPU, there's maybe another 10 GBps squeezable from there, I'll have to see about that too.

Update #3: GPU at 165 GBps (note that this is on a HD 4850 that has 40% of the memory bandwidth of the GTX 285. If the results scale linearly, GTX 285 would get 391 GBps.) CPU at 27 GBps, but the CPU results are wrong so I need to fix them to see if that's a real score.

Update #4: Fixed the optimized CPU version. It's at 22 GBps, but could probably go to 27 GBps with little work. The upper limit for the CPU version is 42 GBps where all reads come from L1. The L2 limit is less than half the L1 limit, so the optimized version is already beating pure L2 bandwidth.

Update #5: According to this article on iXBT Labs, the L2 bandwidth of the HD 4850 is 384 GBps and the L1 bandwidth is 480 GBps, so those would be the limits for the GPU score. If the CPU scores are any indication, 400 GBps lies within the realms of possibility.

Update #6: Wow, a less math-heavy version of the OpenCL kernel absolutely flies on the CPU. 33 GBps. That's equivalent to 8 cycles per cache line. If the Nehalem was doing 10 cycles per cache line, it might go 20% faster with this version and reach 187 GBps. Anyone with an 8-core 2.93 GHz Nehalem want to give it a go (:?

The difference between the GPU version and the CPU version is that the CPU version is doing the dot product with an accumulator vector and doing the horizontal add at the very end of the loop. This is because a float4 dot product is mul + horizontal add + float add on the CPU, and that's going to take 6 cycles or something. On the other hand acc += v0 * v1 executes in one cycle. The GPU does the dot product in a single cycle and the add in another, so it's just 2 cycles. I tried running the CPU version on the GPU but it strokes it the wrong way and results in the GPU hanging (or just running at 0.1 GBps through a benchmark designed for thousand times that). And the ATI Linux drivers don't do the Windows 7 thing of resetting the GPU if it goes unresponsive for several seconds.

There's also a weird data alignment issue (or something) on the GPU: 496x496 images run at 152 GBps, 500x500 images run at 145 GBps, 512x512 runs at 165 GBps. The CPU doesn't care all that much.

Update #7: On trying to optimize the GPU kernel, managed to get the CPU kernel up to 42 GBps by refactoring a bit. That's nuts. That's a Core 2 Duo 2.12 GHz from four years ago beating the paper's two-thread Nehalem 2.93 GHz version by 13%. Maybe I should write a paper titled "Believe it or Not! 2 GHz Core 2 CPUs can Match 3 GHz Nehalems and POWER7s for FLOP-intensive Applications!"

Update #8: Did a few bogus GPU kernels to figure out where the realistic limits are. If the 160 GBps kernel did all its reads from L1, it could achieve 270 GBps.

Update #9: GPU kernel at 182 GBps, achieved by breaking each row into its own extra work unit. Fixed the algorithm as well.. The paper was offsetting both images on Y but the idea is to keep one image static and move the other image over it. Anyhow, that had no effect on performance. It's funny how an extra 15 GBps on the GPU is kinda meh, but a 7 GBps increase on the CPU is WHOAH AWESOMESAUCE. Hurr durr, still lacking the L1 breakthrough. If the new kernel did all its reads from L1, 350 GBps.

Update #10: Oh right, += isn't atomic. Now I get to rewrite things.

Update #11: Actually, += not being atomic wasn't the problem. The problem was that the HD 4850 OpenCL implementation only gives you 16k of local mem. And that's across all the workgroups. And I can't get it to work. Oh well. Anyhow, optimized work group size and got 190 GBps from the GPU (including data transfer to-from main mem), and the kernel bandwidth broke 200 GBps! Woop!

Update #12: CPU kernel at *drumroll* 78 GBps! Insane! That's equivalent to a single quad-core Nehalem socket. From a previous-gen chip with 66% of the clock-speed, half as many cores and quarter of the cache. Applying the same optimizations to the GPU kernel also got a nice 30 GBps performance boost, bringing it up to 235 GBps.

Update #13: Tuned the GPU kernel a bit and got 250 GBPs at 640x640 problem size. Kernel is doing 260 GBps internally. Still need to eke out 16 GBps and it'd match the paper's Power7 numbers. Man, that thing is a beast. Or maybe they were using a version that was more optimized for it than for the Nehalem. My CPU kernel should do, hmm, maybe 500 GBps on the Nehalem, if it scales in a linear fashion. Who knows. [Update: 92% of an 8-core Nehalem's 370 GBps L1 bandwidth would be 340 GBps.]

Update #14: Fixed the kernels, as I had optimized correctness away. 6 GBps drop for the CPU, GPU speed dropped 10-20 GBps.

Update #15: CPU kernel at 84 GBps, GPU kernel at 220 GBps. The CPU kernel is now about as fast as the naive GPU kernel, but still 3x slower than the optimized one. And the CPU kernel produces different results from the GPU kernel and the reference implementation because it does the additions in a different order.

Update #16: GPU kernel at 207 GBps (500x500), 245 GBps (512x512), 274 GBps (640x640), with internal kernel bandwidth reaching 285 GBps on 640x640. The 640x640 bandwidth is about the same as the Power7 number in the paper, yay!

Update #17: 500x500 kernel at 243 GBps. Also, before I was removing the context creation times and kernel build times from the OpenCL timings, but not the context and kernel release times. Which meant that the results weren't really descriptive of real multi-image run performance as both the creation and release are done one time only regardless of the amount of images you process. The creation overhead is 0.4 s and the release overhead is around 0.3 s, which will skew numbers for a test that takes ~1 s to run. To make the times more accurate for multi-image runs, I'm now deducting both the creation and release overheads from the OpenCL runtimes. The runtimes do include the data read and write overheads, as those are image-specific. (Hey, now all my OpenCL implementations have 0.3 s shorter times, yay?)

Update #18: Wrote a GLSL version of the GPU kernel to see if it might be able to use the texture caches better. Short answer: Yes. 204 GBps on a naive kernel.

Update #19: GLSL version at 337 GBps for 500x500 images (350 GBps on 640x640). 70% of L1 bandwidth. Cached memory reads sure are awesome. The performance improvement came from reading in 2x2 blocks. The reason why it's faster is that a texture memory read also loads the adjacent values to cache. It might also be possible to use a similar load-reusing algorithm as I have on the CPU and the OpenCL implementations and make the GLSL version faster. But GLSL is painful for general computation, so maybe not.

The paper

By my calculations^[1] they achieved 160 GBps read bandwidth from the GTX 285 (advertised memory bandwidth 159 GBps), 150 GBps from the Nehalem (advertised peak 2x22 GBps) and 276 GBps from the eight-core Power7 (advertised peak 100 GBps). The algorithm is doing 1 cycle of computation (a SSE mul-add) per 32 bytes of memory read, so it's pretty much a pure bandwidth benchmark. As the total problem size is 8.25 MB, it snugly fits inside the 8.38 MB L3 caches of the CPUs, so maybe the CPU L3 caches are made out of infinite bandwidth magic?

But, uh, the L3 latency on the Nehalem is 39 cycles, which'd give it 37.5 GBps bandwidth at 8 cores. 150 GBps is 10 cycles per cache line. The L2 latency on the Nehalem is 11 cycles. Plus you actually have to do computation for at least four cycles per cache line, though that's going to be masked by prefetch. How can it get 150 GBps?

The paper had a table that compared the different stages of optimization. A single-threaded scalar version ran in 30 s, compared to a single-threaded SSE version that ran in 12 s. A 16-threaded SSE version took 1.8 seconds, which is 6.7x faster than the single-threaded one. Pretty nice with 8 cores. But this is a memory benchmark. How does a memory benchmark get 7x faster by increasing the thread count? That's only going to happen if each core has its own copy of the memory block it's working on. The same thing happened on the Power7 implementation: SIMD + OpenMP resulted in a 14x performance boost(!) So maybe their algorithm is really good at caching, who knows. There's no compilable source code in the paper so...

I did a quick reimplementation of the algorithm in C++ and got 69 s for the scalar version and 58 s for the SSE version on my Core 2 Duo 2.12 GHz. Parallelizing it by tossing a '#pragma omp parallel for' above the top-level loop brought the SSE execution time down to 49 s. Why so slow? The 8 MB input images don't fit in the CPU cache, turning the program from a cache bandwidth benchmark to a memory bandwidth one. By quick calculation, it's plowing through the benchmark at 5.6 GBps.

[Update: A more cache-optimal memory access pattern got the CPU up to 84 GBps, so it's definitely possible to do most reads from L1.]

Now, I also tested the algorithm on uninitialized allocations. And the results were very strange. On uninitialized data, the SSE version completed in 13 seconds. The parallel version ran in 7.5 seconds, for a bandwidth of 37.5 GBps. That's a 7 cycle latency. L2 latency is 14 cycles. It's actually reading from L1.

So how can two 4 MB images fit into 32 kB of L1 cache? The easy answer? They don't. I think what we're seeing here is the Linux VM subsystem doing optimistic allocation. That is, when you allocate a block of memory, all 4 kiB pages in that block are going to map to The Empty Page at first. The pages are only actually created when you write on them. So as all the pages of a 4 MB image actually map to the same empty page of RAM, the benchmark is going to act as if the whole image is in L1.

I next wanted to check out the GPU numbers and wrote an OpenCL version of the algorithm. Which also gave me an excuse to learn OpenCL. I downloaded the ATI Stream SDK and a couple hours later had a working OpenCL version. There's no Linux installer in the package, so you get to do a strange dance setting LD_LIBRARY_PATH and ATISTREAMSDKROOT environment variables.

Writing the kernel was simple. I pretty much copy-pasted the C++ version of the algorithm to do the OpenCL kernel (there are two things that differ between the versions: vec4 turned into float4 and offset loop counters are now job ids). The state setup one has to do to execute the kernel was quite tedious though, and weighed in at 30 lines of code.

The nice thing about OpenCL is that you can also run the kernels on the CPU. And, well, the CPU build of the OpenCL kernel actually ran faster than my parallelized SSE version, which is pretty awesome. The bad thing about OpenCL is that it's like C and OpenGL when it comes to error reporting: it fails silently and you have to milk the errors out. And when you manage to do an infinite loop, the GPU hangs and you have to cold-boot the computer. Additionally, when you're running something time-consuming on the GPU, it freezes compiz for the duration of the computation and you're left twiddling your thumbs. It would be really REALLY nice if the GPU could do pre-emptive multitasking. I've had to Magic-SysRq-boot my computer 4 times today. It's so Windows 95.

But how about the runtimes for the GPU? My first version took 3.5 seconds on a Radeon HD 4850. 14 times as fast as the CPU implementation, twice as fast as the bogus uninitialized version, not to forget having to account for a 0.6 sec kernel compilation and setup overhead. Caching a copy of the built kernel binary to disk and reusing it brought the GPU overhead down to 0.2 s. The CPU OpenCL overhead was 0.2 s with compilation, 0.1 s with binary kernel.

The HD 4850 has a memory bandwidth of 67 GBps, but the benchmark achieves 91 GBps. Did I make an error in my bandwidth calculations? Or might it cache more of the input image than the CPU? Time to write a small memory benchmark... Yeah. 10.3 seconds to do 69 billion 16-byte reads from a 4 MB buffer, 106 GBps. Guess it caches it better than the CPU (which managed 7.6 GBps.)

By spending the same time optimizing the CPU implementation and the GPU implementation, I managed to further drop the realistic case CPU runtime to 12 s (23 GBps) by making the inner loop update the values in blocks. The intuition there being that you're going to have better cache locality if you're updating a bunch of values that use data from the same area. This kind of a thing could make the Nehalem keep the reads mostly in L2, which might well get it close to 150 GBps. [Update: As noted above, you can keep the reads mostly in L1, which might get the dual Nehalem to ~340 GBps, if the Core 2 results scale.]

Applying the same optimization to the GPU version brought its bandwidth up to 165 GBps. That's 10% faster than the $2600 Nehalem duo. On a $100 graphics card. And that's not even close to the 384 GBps L2 bandwidth limit for the HD 4850, so there is probably quite a bit of headroom there. [Update: Well, if you could use the cache with the HD 4xxx OpenCL implementation...]

To wrap up: Benchmarking a cache-optimized algorithm against an algorithm that reads from the main memory is kinda pointless. Or was the goal of the paper to say "Hey, we've got an awesome optimizing compiler that can reorder loops to be very cache-optimal"? Also, I should remember to initialize buffers when benchmarking. And a $100 graphics card trounces $2600 worth of CPU power in this benchmark. Finally, that $100 graphics card sure could use some damn hardware support for pre-emptive multitasking.

I put the sources on GitHub, go check 'em out.

[1] Estimated bandwidth use for the benchmark:


float_sz = 4
out_size = 250 * 250 = 62500

reads_per_out = 2 * (500-0.5*250) * (500-0.5*250) * 4*float_sz
              = 2 * 375*375 * 4*float_sz
              = 4500000
writes_per_out = float_sz

total_reads = out_size * reads_per_out
total_writes = out_size * writes_per_out

total_bw = out_size * (reads_per_out + writes_per_out)
         = 62500 * (4500000 + 4) 
         = 281 GB

Democracy

2010-06-25T15:51:00.005+03:00

You can only call an economy a market economy if at least 60% of the population can fully participate in it on a daily basis and own property.

By this measure, there are several functioning market economies in the world.

You can only call a country democratic if at least 60% of the population can fully participate and vote in a parliamentary meeting about the matters of the state at least once every two weeks.

By this measure, there are zero democracies in the world.

Suppose that you had an economic system similar to a representational republic. Only the members of the parliament are allowed to own property and participate in the economy. Every four years, the people vote to select a couple hundred people who are allowed to participate in the economy. Does that sound sane? A good system of economy?

I blame Socrates and Plato for the travesty that is known as the Republic. They were children of a democracy^[1] and wanted to overthrow the system to install a philosopher dictatorship in its place. The only reason their writings have survived is that they were written in a democracy (and hence not summarily burned) and that they have provided a nice argument for rulers ever since.

"We actually have the best system of governance. Only us well-born/popular/rich/divine/ideologically correct should be the ones to have any political power. See, even the Ancient Greeks who lived in an oh-so-terrible democracy agree: the perfect political system is one governed by a perfect philosopher-king who is omniscient and makes perfect decisions. Not one that's governed by flawed people who know nothing and make imperfect decisions. QED."

Compare:

"We actually have the best system of economy. Only us well-born/popular/rich/divine/ideologically correct should be the ones to own any property. See, even the Ancient Greeks who lived in an oh-so-terrible free market agree: the perfect economic system is one governed by a perfect economist-king who is omniscient and makes perfect trades. Not one where property can be owned by flawed people who know nothing and make imperfect trades. QED."

[1] Athens had a population of around 200 000, of which only 30 000 were members of parliament, so it fails to reach the 60% standard. However, that was 2500 years ago, in the times of the Persian and Egyptian God-Kings, the Spartan military dictatorship, the Warring States Period, etc. and should perhaps be judged in that context. In a modern society, where you have equality of sexes and no slavery, they should have had 120 000 members of parliament. Then again, where are the modern countries with even a paltry 30 000 members of parliament.

Funny multiplication tricks

2010-06-24T02:11:00.005+03:00

This kinda struck me yesterday: multiplication is additive pattern splatting (I just came up with that term.) I've been seeing binary tricks that have used multiplication as copy operation but never really understood it before now.

With sparse binary patterns it's quite apparent. For example, let's multiply 00100010001001001 by 111. Sounds hard? It's actually easy, you replace each 001 in the long pattern with the short pattern and end up with 11101110111111111. In the same way you can copy a byte to different places in a word. This may sometimes be faster and/or less bother than bit-shifting:


int64 fillInt64WithByte_shift(int8 b) {
  return b << 56 | b << 48 | b << 40 | b << 32 | b << 24 | b << 16 | b << 8 | b;
}

vs.

int64 fillInt64WithByte_mul(int8 b) {
  return 0x0101010101010101 * b;
}

It also works with decimal numbers: 10101 x 99 is 999999. If you have a number that's not a 1 in the target pattern, you need to multiply your source pattern with that: 001 002 003 0005 x 203 = 203 406 609 1015.

However, the simple copy splatting only works when you don't have overlaps. With overlaps you need to add the new splat to the target pattern. Let's consider non-carry addition first using 1001301001 x 211 as the example. First you fill in the non-overlapping parts (the 1001001001 pattern) and end up with 211 211 211 211. To integrate the 3 into the result you multiply it by 211 and add the result to the target pattern.


   1001301001 x 211 =
 211211211211
+    633
 ------------
 211274511211

When you have carries, it gets a good deal more tedious. Consider binary 1111 x 1111. First we splat 1111 on the first 1 in the target pattern: 01111000, then on the second 1: 00111100, etc. and finally get 01111000 + 00111100 + 00011110 + 00001111 = 11100001. Drudgery.

We could move it to base-15 and get 10 x 10 = 100, then convert that back to base-2... for no benefit whatsoever. 1111^2 x 1 = 1111 x 1111. Or do a trick with 10000 x 1111 - 1111 = 11110000 - 1111 = 11100001.

Or convert to base-16 for F x F = E1, then convert that to base-2 for 1110 0001. Note the similarity between F x F = E1, 9 x 9 = 81, o7 x o7 = o61, b1 x b1 = b01. This can be explained by


For base N, N = 10 - 1:
N x N =  N x (10 - 1)
      =  N x 10 - N
      =  N x 10 - (10 - 1)
      = (N x 10 - 10) + 1
      = (N-1) x 10 + 1

Back to splat-additive multiplication, it seems to work for positive integers. Does it work for negative ones? Well, yes, 1001001 * -101 = -101000000 + -101000 + -101 = -101101101. And there's that sign-trick to turn multiplication of negative integers to multiplication of positive integers: -A x B = -(A x B), -A x -B = -(-(AxB)) = AxB.

Does it work for matrices... not that I can tell. Matrix multiplication is not splatting anyhow. What would a splat-additive matrix multiplication operation look like? I don't know, let's start exploring by trying to extend the scalar splat-additive multiplication to vectors. If we treat vectors as bucket-based integers (that is, no carries, base-∞), integer-like splat-additive vector multiplication would be quite easy. For example, [1 0 0 2] [x] [3 6] = [3 6 0 6 12], but what can that be used for? Is there some sensible conversion to scalar from it? Can it be extended to matrices?

As an aside, scalar vector multiplication can be used to make integer multiplication a bit simpler. Using 32 x 1935 as an example: 32 x [1 9 3 5] = [32 288 96 160]. We can convert that to an integer value by doing a base conversion for each element and adding them together: 32*1000 + 288*100 + 96*10 + 160*1 = 61920 = 32 x 1935.

Another way to do splat multiplication with vectors would be to generate a matrix from vector multiplication, an approach not entirely without appeal. [A B] x [C D] would result in [ [AC AD] [BC BD] ]. It can express number multiplication as a special case: AxB would be |A| x |B| = AxB (I'm using |A| to say that |A| is a single-element vector of zero dimension). If you have dimensional numbers, it preserves the dimensionality of the result, [Length] x [Width] = [[Area]].

Can the vector-splatting multiplication be extended to vectors of vectors? Seems so:


[ [1 2] [3] ] x [ [5 5] [6 2] ] would be, uh, ... 
[  [   [1 2]x[5 5]     [1 2]x[6 2]   ]  [   [3]x[5 5]  [3]x[6 2] ]  ] =
[  [  [[5 5] [10 10]] [[6 2] [12 4]] ]  [ [[15] [15]] [[18] [6]] ]  ]

And to continue with the dimensional numbers,


[Height] x [[Area]] = [[[ Height x Area ]]]
[[Area]] x [[Area]] = [[[[ Hypercube! ]]]]

What properties would this vector-splat multiplication have?


[ A B ] x [ C D ] = [ [AC AD] [BC BD] ]
[ C D ] x [ A B ] = [ [CA CB] [DA DB] ]
Not commutative (unless the equality is sum-equality: A = B if sum(A) = sum(B).)

([ A B ] x [ C D ]) x [ E F ] = 
[ [AC AD] [BC BD] ] x [ E F ] =
[ [[ACE ACF] [ADE ADF]] [[BCE BCF] [BDE BDF]] ]

[ A B ] x ([ C D ] x [ E F ]) =
[ A B ] x [ [CE CF] [DE DF] ] =
[ [[ACE ACF] [ADE ADF]] [[BCE BCF] [BDE BDF]] ]
Seems to be associative.

The identity element would be a dimensionless 1:
1 x [ A B ] = [ 1xA 1xB ] = [ A B ]

Note that a dimensional one won't do:
[1] x [ A B ] = [ [A] [B] ]

To do the inverse element, we need vectors with negative dimensions:
I x [ A B ] = 1 requires I to reduce the dimension of [ A B ]

Additionally, the inverse element requires sum-equality (i.e. [A B] = [B A])
[A B] = [C] iff A + B = C
|0.5 0.5| = 0.5 + 0.5 = 1, hello confusing notation for dimensionless vectors
[[A [B] C]] = |[[[B]]] [[A C]]|, hello polynomial notation

I x [ A B ] = 1
[1/2A 1/2B]^-1 x [ A B ] = |A/2A B/2B| = |1/2 1/2| = 1

So the inverse function would be
-- get inverse vector for which V x inv(V) = 1
inv(V) = (map (\i -> 1/(V_len*i)) V)^-V_dim

...ghhh. So, uh, I guess with sum-equality it could maybe form a field with some sort of vector addition. The real question is whether there is anything where this multidimensional splat-multiplication works better than existing methods. It looks a lot like polynomials. Maybe the vector addition could be..


|A B| = (A + B)₀
[A B] = (A + B)₁
|[[A B] C] D| = |[[A B]] [C] D| = (A+B)₂ + C₁ + D₀
= [[A B]] + [C] + |D|

That's kinda interesting, the (A+B)₂ isn't (A+B)² but rather means that (A+B) is a two-dimensional number.

Oh well, back to the grind.

Status update

2010-06-23T01:19:00.003+03:00

Taking two days off, a mid-week weekend.

Estimated that doing an editable zoomable tilemap would take a week. 15 days of consecutive commits say that I was wrong by at least a factor of three. Zoomable tilemap, check. Editable, sort of. Lacking in genericity as rendering in browsers is either slow or requires more ugly performance hacks than you can shake a stick at. And 60 fps is the law of the land, so screw genericity.

OCaml is horrible on 32-bit. 16MB maximum string size. And strings are the intended way to handle binary data. And if you use BigArrays, they're not GC'd properly. Then again, the worst-case allocation in my use would be 2.5 gigs, so perhaps it would be reasonable to work very hard to do things in smaller increments. Loading 16k x 16k images, that is my nemesis (I had a Ruby + C version of this but...)

Next need to either make an interactive character animation system in JS or find one. Hrm.. Maybe look at the Big Buck Bunny Blender files, see how to do the whole model-rig-animate-thing, export to JS, write some simple engine to run it. Or just draw the animation by hand (On paper! With a reed! A broken reed!)

WebGL stereo rendering demo

2010-06-06T17:33:00.002+03:00

Cross your eyes and prepare for eyestrain, here comes the WebGL stereogram. Drag with LMB to rotate, drag with MMB to pan, wheel zooms, shift-wheel scales model.

The demo uses glScissor and glViewport to restrict drawing to one side of the viewport, then draws the scene to the left and right sides with the camera offset slightly. The camera offset is along the right-vector of the camera and you can get it by finding the normal of the camera up and look vectors. Basically right = cross((lookAt - position), up), then leftEye = position - (sep/2 * right); rightEye = position + (sep/2 * right). Then you render to left viewport from rightEye and to the right viewport from leftEye.

Loading WebGL models from tarballs

2010-06-02T20:38:00.005+03:00

Plugged the TarGZ loader together with my old OBJ and binary model loaders and emerged with a couple demos that load a model and a texture from a gzipped tarball and render it using WebGL.

Verdict: use uncompressed tarballs and HTTP gzip compression whenever possible. Doing the gzip with JS doesn't cache the uncompressed data transparently (you'd have to do an offline storage hack), and it's kinda slow. Parsing tar vs. plain files is a bit different tradeoff. With tar you only have to manage a single file. With plain files, you need to manage the whole directory hierarchy they live in. Which may not be something that you can or want to do with your file serving host.

One solution is uncompressed tarballs and a .htaccess with deflate filter set for them, though it seems to screw up caching somehow (thanks for the .htaccess tip go to Evgeny Demidov)


AddOutputFilterByType DEFLATE text/html text/plain application/x-tar application/javascript

<Files *.bin>
SetOutputFilter DEFLATE
</Files>

which also helps for JSON models and OBJs and other textual things. The main problem with the tarballs is that loading the images takes a bit longer (read: hundredths / tenths of a second) due to having to dataURLize them. Well. If you pre-dataURL the images, that won't be a problem.

The main benefit of the tarballs is that everything's in a single file and you don't need a chain loader to initialize models in a single step. Instead of first showing an untextured model, then sometime later the model with a diffuse texture, then later the model with a diffuse texture and a normal map, etc., with single-step initialization you get the full textured model once all the data is loaded. Plus, if you have a standard model format of some sort with manifests and whatnot, you can swap models by changing a single URL. You can do that with a plain model.xml or somesuch as well, but then you need to manage multiple files. It's harder to reuse resources inside tarballs than resources as separate files, though.

Also updated the big whale model viewer to be mouse-controllable (drag with LMB to rotate, drag with MMB to pan, scroll wheel to zoom, shift-scroll wheel to scale).

Holographic video

2010-05-31T18:45:00.007+03:00

NHK is developing a holographic TV that doesn't require viewing glasses (via Reddit.)

The premise is simple: put a transparent surface between the viewer and the object. All the viewer sees is the photons coming from the transparent surface. Now sustain the photon emission from the transparent surface as-is and remove the object. To the viewer it looks like nothing has changed. The object seems to be still visible through the surface and you can walk around and look at it from different angles. A hologram! (When I understood the idea, I got a stupid goofy grin and walked around the house going all "WOOHOO! HOLOGRAMS! YEAH!")

Implementation then. The camera lens could be a large sheet of viewing-angle-optimized hemispherical pixels with an RGB light detector for each direction. The viewing device would swap the detectors for emitters. It'd require a huge lens (the size of the display, really) and require multiple times the megapixels of the usual one-directional image (as now each pixel is in fact several pixels, each facing a different direction.)

Except that they're not doing it quite like that. I don't fully understand it, but I think they're using an 2D array of lenses to get a flat projection and then record the resulting 2D image in the usual manner (project to a small CCD/CMOS/whatever). The display is a regular ultra-high-resolution flat panel with a similar array of lenses in front of it. Here's a good presentation slide deck: 3D-TV System Based on Spatial Imaging Method for Future.

Loading .tar.gz with JavaScript

2010-05-28T20:52:00.002+03:00

With the power of copy-pasting code from the INTERNET, I have now created a JavaScript library that streams gzip files and specifically .tar.gz files.

DEFLATE and CRC32 implementations found with Google and tweaked and bolted onto the tar parser. The performance is something between 1-5 MB/s on Chrome, third of that on Minefield.

In action: Demo page
Just the library: gzip.js

Numbers, numbers

2010-05-25T15:57:00.004+03:00

A Finnish copyright survey tells us that copyright infringement costs the Finnish recording industry more than 355 million euros every year. The yearly sales of the Finnish recording industry are around 100 million euros. If 30% of that is profit — 30% being the profit margin for software companies — they spend 70 million to make the records, get 100 million back in sales and pocket a cool 30 million in profit.

Now, as a result of the copyright infringement costs, the costs of the Finnish recording industry are actually 70 + 355 = 425 million euros. They spend 425 million euros, get 100 million back in sales and end up 325 million in the red. I will expect the whole industry to go bankrupt in three months, if not faster. Every euro they spend is costing them five.

Slideshow-qgl, now with new features

2010-05-24T19:12:00.006+03:00

Slideshow-qgl — my tiny slideshow app — now has recursive dir traversal, a fullscreen mode and a slideshow mode. Which makes it more useful than.. most slideshow apps out there. GQview and EOG, I'm looking at you. With a stern eye. And a hint of rebuke.

As you might guess from the name, slideshow-qgl is based on Qt4 and OpenGL. Qt4 makes it relatively easy to write and OpenGL makes it fast to pan & zoom & fade between images.

It follows the mplayer approach to GUI design, specifically "widgets are ugly, get rid of them." This is what it looks like with an image loaded and zoomed in a bit:

Right.

Ugly Tachikomaling

2010-05-23T02:22:00.000+03:00

Parsing tarballs with JavaScript

2010-05-21T09:09:00.024+03:00

Update: Check out this augmented version that streams gzipped tarballs.

Here's a small piece of JavaScript to parse tarballs and my custom JSON packfiles. There be four demos as well: loading files from a tar, streaming images from a tar, loading files from a JSON packfile and streaming files from a JSON packfile.

The part that converts images to date URLs is a bit slower than it could be, as it has to strip high bytes off the characters. The upcoming JS File and Blob APIs for binary data handling should help there. Though if you have less than a hundred kB of images, I don't think you'll even notice the delay. Even half a meg of stuff unpacks in a fraction of second on my slow laptop (Pentium M 1.7GHz). If you do need speed, you can convert the images to data URIs beforehand.

Quickly estimating, it'd take something like fifteen seconds to load up a hundred megs of models and textures on my laptop, maybe around 5 s on a decent computer. Doing the initial archive parsing pass would take maybe a second for a hundred meg archive. If that's too slow for you, I want your internet connection. If the hundred meg tarball is split 1:4 geometry:textures, where the geometry takes 20 bytes per tri and the textures are 10x compressed JPEGs, it'd have 1 Mtri geometry and 240 Mpx textures.

The script doesn't handle gzip or any other compression, use gzip-encoding on the server for that. The tar file format is pretty simple: it's based on 512-byte blocks and each file begins with a 512-byte header, followed by the file data padded up to a multiple of 512 bytes. The numbers are represented as octal ASCII (though there is a GNU tar extension that uses binary ints for handling files bigger than 8 GB, which my script doesn't support).

My JSON packfile format consists of a one-line JSON header array of {filename : string, offset : bytes, length : bytes} followed by a newline and the concatenated file contents. Easy to create and parse.

Edit: Added streaming using xhr.readyState == 3 checks. It might cause some stuttering on the page when dataURLing the images, though it should be quite efficient otherwise. Optimizations welcome :)

Vanishing hardware

2010-05-20T23:58:00.006+03:00

... Apparently the ethernet chip on my desktop computer motherboard completely disintegrated and now it doesn't even show up on lspci. All I know is that it disappeared after waking up from Windows 7 sleep mode.

Oh well, another excuse to upgrade.

Edit: after putting Windows to sleep again and waking it up, it seems that it momentarily lost and then regained the SMBus controller as well (not the ethernet controller however), in the process deactivating Windows. How very Windows of it. A pox on it! A pox on it and its descendants to the seventh generation!

Energy use by economic output

2010-05-20T01:35:00.003+03:00

Energy use in kWh per GDP per capita. Country data from NationMaster, plotted with R. The green line is the LOESS regression curve. The blue lines are averages for top and bottom 100 countries (after removing the best and worst 10 countries), the pink lines are averages for countries ranked 11-20 from top and bottom.

Roughly speaking, the general trend is five bucks GDP for every kWh of energy use. The least efficient countries get 2 $/kWh and the most efficient 25 $/kWh.

Assuming that energy use follows economic growth and that the average annual global economic growth is 3%, we can extrapolate (hah) that the global energy use would equal the total output of the Sun at around the year 3125. Put that in your pipe and smoke it.

A chain of trivia

2010-05-17T18:48:00.010+03:00

Chemical makeup of the human body - Elements by mass 65% O, 18% C, 10% H, 3% N, 4% others. You're three times as oxygen-rich as the air you breathe. 98.7% of the molecules in your cells are water molecules (but only 65% of mass, water molecules are light). There are around 174 teramolecules of water per cell. Your DNA chromosomes are huge giant molecules about 5 billion times more massive than water molecules.

R-generated graphs of the first fifteen elements in the human body by mass and count. Hey, power law distributions!

Elements in Animals and Humans - what are the other elements used for? Deficiency in Magnesium causes tetany, muscle spasms.

Clostridium tetani - the Tetanus is an anaerobic soil bacteria that seems to produce one of the most deadly neurotoxins for no good reason whatsoever. The Tetanus toxin binds to muscle nerve synapses and blocks the release of relaxation signals, causing death by bone-shattering spasms.

C. tetani is closely related to Clostridium botulinum, another lethal anaerobic soil bacteria. The Botulinum toxin (a.k.a. Botox) binds to muscle nerve synapses and blocks the release of activation signals, causing death by flaccid paralysis.

The two above bacteria are Gram-positive, as opposed to Gram-negative, the difference between the two types being that Gram-positive bacteria retain a violet crystal stain whereas the Gram-negative ones can be re-stained with a different color. The cause of this is apparently that Gram-negative bacteria have an outer layer around their cell walls that doesn't absorb the violet crystal stain as readily.

Animal cells and protozoans don't have cell walls. Plant cell walls are made out of cellulose, bacteria cell walls consist of peptidoglycan. The antibiotic enzymes, Lysozymes, react with peptidoglycan, breaking a molecular bond in it to tear apart the bacterial cell wall (by contorting a sugar in it to an uncomfortable position that breaks easily). Fungi cell walls are chitin, as are the exoskeletons of insects. Diatoms, those strange silicon phytoplankton of the seas, have cell walls made out of silicic acid. Diatom shells along with Radiolarian shells form the siliceous ooze at the bottom of the sea.

Pelagic sediments are the result of the constant rain of sediment over the abyssal plains. The rain consists of dead plankton shells, remains of marine animals, micrometeorites and dust carried by rivers and winds. If you ever wondered where marble comes from, the chain starts here. Calcareous plankton shells fall to the ocean floor, forming calcareous ooze, the precursor of limestone, a type of sedimentary rock that consists of calcite (CaCO₃). The limestone then gets hammered in the forge of crustal plate collisions and recrystallizes to marble, some of which was dug up in Attica at around 440 BCE and used to build the Parthenon, the CaCO₃ of which is now being dissolved to aqueous Ca²⁺, SO₄^2-, H₂O and CO₂ in a reaction with sulfuric acid (H₂SO₄) from acid rain.

The clouds covering Venus are composed of sulfuric acid and it rains sulfuric acid in the upper atmosphere, but the rain evaporates around 25 km above the surface where the temperature hits 300 °C. The surface temperature of Venus is around 470 °C, which is close to Draper point (525 °C) where solid materials emit clearly visible blackbody radiation. Maybe dark places on Venus could glow a faint dull red.

The deep-sea hydrothermal vents glow as well, but with a spectrum too blue to explain with just blackbody radiation from the hot 350 °C water. The higher frequencies might be from perturbation caused by the dynamic environment, namely sonoluminescence from tiny bubbles bursting and crystallo- and triboluminescence from crystals forming and breaking.

Athenian Democracy

2010-05-16T04:41:00.008+03:00

Listening to this Athenian Democracy lecture by Donald Kagan. The Wikipedia article is detailed too.

Ancient Athens of the time had 100,000 citizens, of whom maybe 30,000 were politically eligible (male citizens of 20+ years of age). They had political meetings every ten days to discuss and set policy about everything with majority votes. These political meetings had maybe 5,000 people present each time, but there were some votes that required 6,000 or more present to be binding.

Compare that to a modern representational democracy such as Finland. 6 million citizens, around 4.2 million eligible (citizens of 18+ years of age). They vote once every four years to choose which parties should get parliamentary seats. The parliament numbers 200 and meets four times a week to set policy about everything with majority votes.

If the Athenians had used a similar system (scaled by population), the number of members in their parliament would have been one.

Quick game reviews

2010-05-15T03:58:00.006+03:00

The night is dark and long and I have nothing to do.

Puzzle Quest, a character-building story battle game with movement done on a map and battles resolved in Bejeweled. 65 hours played. Addictive, not all that fun to play, the match-three matches take too long, the item & spell combos are fun to come up with but haven't gotten powerful enough to shorten the matches. The art is nice but the gameplay GUI art is not so nice. Story is whatever but has a non-stereotypical stereotypical fantasy world with the player helping everyone regardless of what they want to do. And what they want to do always seems to involve playing Bejeweled against megafauna / Bob Warlord. Playable in 15-min increments in the same way that you can smoke just one cig a day. Story: A bad guy is resuscitating people and no-one wants to help and minotaurs are turning into robots to kidnap dragons and orcs and ogres are being belligerent and you must put some other minotaur bad guy together and resuscitate it so that it can fight the other bad guy and my head is hurting.

Company of Heroes, a WW2 RTS about the U.S. experience on the Western Front to take out of the remaining 20% of the Nazi army. 19 hours played. Nice graphics, jaggedy shadow maps, paintings shown between missions lack finish. Battles take too long, three guys building a tank factory to build an endless supply of Shermans makes it very RTS gamey. Not casual, takes hours per campaign map. Cover system good. UI less good than in SupCom. Not all that fun to play, but not addictive either. Story: terrorists have taken over Europe and you must kill them all. I guess?

The Misadventures of P.B. Winterbottom, a cloning-based puzzle platformer about pies. One hour played. Good presentation. Strong style. Good music loops. Figuring out the tricks to the levels is fun but means frustration, which is less fun. Not quite casual, requires ~half hour per session. Story very wtf.

Penny Arcade Adventures Episode 2, a character-building story battle game with movement done by walking characters around sets and battles resolved in turn-based combat. Half an hour played (plus Episode 1). Much like Episode 1. Maybe too much. Gabe's special attack minigame is no longer the most frustrating thing ever, which is good. Good presentation, strong style. Has battle-helping items with a time penalty to utilize, making one less likely to bother. Not quite casual, half hour sessions. Story: there are more robots and you should break them all.

Star Wars: Knights of the Old Republic, a character-building story battle game with movement done by walking characters around sets and battles resolved in pausable real-time combat. 1.5 hours played. Good presentation. Oh god do I have to talk to everyone here and beat up military people and this racist city is very boring why do I have a sword and a glow-bolt pistol where should I go now this is a lot like Baldur's Gate I'm assigning skill points to some skills I don't know anything about, hey I'm dying in these combat sequences, how quaint. I'll play more at a later date (said five months ago.) I guess it has a story to it?

Torchlight, a character-building story battle game with movement done by walking characters around sets and battles resolved in real-time combat. 41 hours played. Nice presentation, good ambient music, randomly generated maps, way too many drops, decent feel of POWER in attacks (I would still take it up a notch but I'm weird like that), fiddly targeting. Addictive. If Diablo 2 is Zangband (well, noooot really but), this is Moria. Go down down a dungeon, down down down. Down. When done with the story dungeon, there's a new random dungeon that never ends. The random dungeon generation system is nifty. Summoning golems is nifty. Lightnings and laser beams are nifty. Can be played in short stints but takes a while to load and is kinda iffy windowed. Story: There's a McGuffin in the earth and it turns you into a monster. Except not really. Only the bad guy can do that. In the meantime you just genocide everything living down there because they're hindering your walk.

Warhammer 40k: Dawn of War 2, a character-building story battle game with movement done by clicking on a map and battles resolved in tactical RTS-ish real-time combat. Played 20 hours. Randomly generated maps (I.. think.) Boss fight at the end of each map. Playing it feels like using Photoshop, very much a keyboard shortcut game. Move your four squads of bald robot men from cover to cover in a very micro fashion. Starts off challenging, gets difficult, mocks you for doing less than perfectly. Game gets uglier the more you screw up (invading monster army turns the air green with red flakes). Not casual, maps take hours to play and you're punished for failing. Presentation much like Company of Heroes (same engine, same studio) but with more shiny metal bits and glowy glow guns. Story: an invasive space species carried by the space winds is taking over space planets and the indigenous surface dwellers try to chase it off before it turns their C-H-O into its C-H-O.

CoH vs. DoW2: In WW2 you try to look like a tree or a grass because you are a hunter and don't want to be spotted. In the grim future where there is only war you wear colorful shiny clothes because war is waged as an amusement and proving ground for the aristocrats, so you must have bright clothes to show off the colors of your team.

Trees and grass are a form of water and air that absorb solar radiation to turn more water and air into tree and grass. Hence they try very hard to not be shiny reflective, hence clothes that try to look like they're trees and grass are also not shiny, hence a WW2 infantry game is not shiny. Also, you're more likely capable of digesting plastic than metal. Because plastic is made of C-H-O, same as you. Metal is not what you're made out of, so you can't turn it into your own mass, so you don't have a way to digest it. And you're mostly made out of air and water too. With some small gunks of dissolved crust elements such as Ca, P, K, S, Cl, Na, Mg, Fe.

Beyond Good & Evil, a story game with photography and occasional real-time combat, movement by walking characters on sets. Two hours played. Good presentation. Strong style. Controls on PC quite horrible. Performance problems. Tells a story, creates a world. Should finish it.

Trine, a character-building physics puzzle platformer with real-time combat. 1.5 hours played. Very nice level graphics and main character models. Paintings less nice. Somewhat frustrating. Story whatever at least up to now. Puzzle Quest ate all my time so I haven't played this much. Story: A McGuffin spares game designer from having to come up with a way to have three separate characters working together in a single-player action game.

Prince of Persia (2008), a running and jumping game with real-time combat. Played it through. Production values through the roof. Fun to play, pretty scenery, good sounds, strong style, non-obtrusive storytelling. You're constantly being saved. Story: bad gods help people, good gods do nothing, felling trees is bad but sometimes you gotta do it.

Hammerfight, a skill-based real-time combat game. Played the demo through. Super difficult, but good impact when your attacks connect. GUI graphics rather awesome, game graphics good too with glowing colors and banners waving in the wind. Slow renderer, most lacking in fillrate for smoke. Talking heads in story not so good. Story: Neighbouring emperor wants to conquer you! And he succeeds! Man.

Why do I buy these games if I don't play them? Probably because it's easy and they're cheap. And grocery shopping desensitizes you to paying 10e for something.

Git commit stats

2010-05-13T23:48:00.008+03:00

Tonight I wrote a little tool in Ruby to print out stats based on git log output. See the project page on GitHub. (Hey, apparently there's a ready-made tool for this. And it has the same name too! Better change the name of my script to gitnuts.)

Here are the weird things it tells about my commits on a project:


                  Commits: 1365
            Changed files: 2483
 Changed files per commit: 1.8
        Total added lines: 86277
      Total deleted lines: 22628
      Total changed lines: 108905
         Added per commit: 63.2
       Deleted per commit: 16.6
       Changed per commit: 79.8

Changed lines per commit
Bucket size: 10
Bar step size: 9
      0 *********************************************************** 538
     10 ******************** 187
     20 ***************** 161
     30 ********** 98
     40 ******* 68
     50 ****** 62
     60 **** 43
     70 *** 33
     80 * 14
     90 ** 18
    100 * 17
    110 * 15
    120 ** 21

Commits per month
Bar step size: 7
2008-01 ***** 35
2008-04 **************** 117
2008-05 **************************************************** 367
2008-06 ********************************************************** 411
2008-07 ******************************************** 313
2008-08 **************************************** 286
2008-09 ********************* 148
2008-10 ** 15
2009-01 *** 22
2009-02 **** 32
2009-10 ************* 92
2009-11 *********** 77
2009-12 * 9
2010-01 ****************************************** 295
2010-02 ********** 71
2010-03 ************ 89
2010-04 * 13
2010-05 ********** 76

Commits per weekday
Bar step size: 7
  0-Sun ***************************************** 291
  1-Mon ****************************************** 298
  2-Tue ********************************************* 317
  3-Wed ********************************************************** 406
  4-Thu ********************************************************** 408
  5-Fri ********************************************************* 404
  6-Sat *************************************************** 359

Commits per hour of day
Bar step size: 3
      0 ************************************************ 145
      1 ******************************************************* 165
      2 ************************************ 108
      3 ********************************** 103
      4 ***************************************** 124
      5 **************************************** 121
      6 ************** 42
      7 ****************** 54
      8 ******************* 57
      9 ************ 37
     10 *************** 46
     11 ******************** 62
     12 *************************** 82
     13 ********************** 68
     14 ********************************** 104
     15 *********************************** 107
     16 ************************************ 109
     17 ********************************** 103
     18 *************************************** 118
     19 ********************************************* 135
     20 **************************************** 122
     21 ****************************************** 127
     22 ****************************************************** 164
     23 ************************************************************ 180

How about some chartjunk courtesy of R? Check color identifies commit author, y-axis is log10 of change size + 1, x-axis is time in seconds from epoch. The lines are lowess smoothings of the raw data.

Effort estimation game first round results

2010-05-08T03:44:00.001+03:00

First round results

Task 1 estimated: 50 lines, 30 min implementation, +15 min making it work right.
Task 1 measured: 55 lines, 60 min implementation, +15 min making it work right.
Task 1 adjustment factors: 1.1x lines, 2.0x implementation time, 1.0x fix time.

Task 2 estimated: 10 lines, 10 min implementation, +15 min fix time.
Task 2 measured: 39 lines, 15 min implementation, +12 min fix time.
Task 2 factors: 3.9x lines, 1.5x impl, 0.8x fix time.

Second round

Task 3 estimated: 150 lines, 8 hours implementation, 17 hours fix time.

Task 3 estimates adjusted with factors from tasks 1 & 2:
375 lines, 14 hours implementation, 15 hours fix time.
Let's see how it goes...

Bayesian effort estimation game

2010-05-05T23:59:00.004+03:00

Update: First round results!

How accurately can you estimate how much effort it takes you to write a new piece of code? Fear not, I have a solution! An untested solution! Very academic, you might say!

Compile a list of features of different sizes that you're going to implement.

Pick a small feature for starters.

Estimate the effort required for it: lines of code, time to "it compiles!", time spent fixing bugs, time spent testing, time spent documenting.

Document the effort: Write in your notebook when you start and stop working on the feature, and the lines of code (look at the commit data). If you don't work on the feature, write why.

When you're done with an estimated task, compare the estimate and the effort by computing the effort/estimate -ratio for the task. This is your new a priori estimate error for that task.

Pick a new feature, estimate, multiply the estimate by the a priori estimate error.

When completing the new feature, calculate its effort/estimate -ratio. Add the ratio to the list of ratios, use their average as the new a priori estimate error. When you have a bunch of errors, plug them into R and start keeping track of their mode, standard deviation and whatever useful statistical variables you can think of. Plot a curve for estimate error per task size.

The goal of the game is to keep the average estimate error of the previous 10 tasks as close to 1 as possible as long as possible.

Take photographs of the people playing this game and make trading cards of them with seasonal averages and all the interesting statistics you collected earlier.

Sell the trading cards to kids.

Profit!

Drawing with writing

2010-04-30T20:04:00.009+03:00

The four faces on the left are made up of Hanzi/Kanji radicals, the middle face is hiragana, the left-side faces are Hangeul, Latin, and Greek. The face left of the bunny and the bunny itself are Cyrillic.

One thing that's been interesting me during the writing system study is how writing systems affect the drawings of the people who use them. I mean, if you're going to draw something, you're going to use the strokes and patterns you already know. And the patterns you know best are the ones you use every day, so the writing system should have a impact on the drawings. There seems to be a feedback effect at work here: the things you draw turn into writing and the writing guides your drawing.

To explain, drawing and writing are both ways to communicate information visually. The old writing systems were in fact formalized alphabets of simplified drawings (or not so simplified drawings in the case of Egyptian hieroglyphs). And these drawings were based on the surrounding nature of the time the writing system was made, hence you have cobras, pintail ducks, flamingos, jackals, crocodiles and hippos in the hieroglyphs. Had it been a Northern European writing system, it would've had wolves, elks and bears instead of the Nile fauna.

One example of writing systems affecting drawing that's easy to reason about are emoticons. You mostly do emoticons that are easy to type and don't get mangled in the transmission. With a Latin keyboard layout you get emoticons with Latin characters and punctuation. With a Korean keyboard, you get Korean emoticons. There's a feedback effect here as well: emoticons create an alphabet of emotions, that you then use in your drawings as a shorthand for the emotions in question. And when you draw an emotion without an emoticon, someone might figure out a way to simplify that drawing into an emoticon, which is then added to the alphabet and used by people drawing emotions...

And then you have writing systems designed for drawing. Perhaps I could segue from drawing writing systems to learning circuit analysis.

To take this drawing with writing thing further, Islamic calligraphy is one famous example:

Image courtesy of Wikipedia.

And the Egyptian Book of the Dead... those big guys there are more like 180pt calligraphy than drawings. Even their sculptures veer towards being big 3D letters carved from stone.

Images courtesy of Wikipedia. Click to see the pages.

That hand-wavy movie GUI

2010-04-30T02:00:00.003+03:00

Oblong's G-speak prototype (G stands for Gesture, as far as I can tell). They consulted on the Minority Report movie, which then inspired them to build an actual prototype. (Via Cocktail Party Physics)

Who do you sell that to? Museums? Universities? Put it in lecture halls, have the professors wave through their slide decks? Automotive design houses might buy a couple to augment their VR caves. Air traffic controllers... kind of, except that it needs to be either secondary or a proven stable system.

It is the new technology dilemma: no current applications, new hardware, no exact target market. (Or very entrenched "target markets" that want proven products which have jumped through all the regulatory hoops and are backwards compatible and have an excellent sales force to wine and dine the people holding the strings of the corporate purse.)

I suppose the largest market would be semi-wealthy individuals. Strap it onto a 50" HDTV, have them wow visitors with awesome gesture-driven slideshows of holiday pics and the way you can change channels and surf the web from your couch. But without a remote or a keyboard. Worked for the iPad.

You can't replace the computer desktop with it, because the desktop has infinite inertia. And you need to rewrite all the apps to fit that thing. It does sound like they do need to use a carbon copy of the iPad strategy: have an App Store and integrate media stores and the web, require creation of a Store account to use it, have it be a not-computer (to not set high expectations and to avoid having to compete with the strong competitors), make it look cool and different (again to avoid comparisons), do the cool things well and basic things passably, sell an upgrade a year later ("It fixes all your problems!")

Extract money from people entering the platform, extract money from people on the platform, extract money from people upgrading the platform, extract money from people who want to sell things to people on the platform. Fight very very hard to maintain your grip on the money pipes.

Oh, also, I drew a new picture.

Early (and recent) computer art

2010-04-28T14:43:00.007+03:00

VISUAL AESTHETICS IN EARLY COMPUTING (1950-80) (via Trivium)

It's the demoscene, 50s-style!

Or maybe something akin to modern dance with a webcam-operated CFD backdrop.

Speaking of the demoscene, here are a couple recent productions:

Fairlight Cncd - Agenda Circling Forth. Millions of particles. Nice forest scene.

Farbrausch - fr-043: rove. Good ambience.

United Force & Digital Dynamite - Wir sind Einstein. See also The Golden Path for more trippiness.

I do wonder if we'll ever have user interfaces with those kinds of production values. The recent trend has been to lock the GUI down so hard that you can't even change the theme, sometimes you can't even change the colors (OS X, I'm looking at you.) On the other hand, web browsers seem to be picking up some easy themeability with those background-image themes of Chrome and Firefox. Squeeze creativity out from one place and it goes to another.

But still. There aren't many GUIs that look nice and fun. Well, except in games. In games they're the selling point. The fancier and better your game GUI looks (i.e. the thing we call "graphics"), the better your game will sell. It's all about the experience, the fancy wrapping around the game logic crack cocaine. Sure, you want to _play_ the game, but having the process of doing that look and feel awesome makes the experience so much better.