John Myles White: Die Sudelbücher

In light of human foibles, I would suggest that this decision rule be used in lieu of Occam’s Razor: of several possible explanations for an observation, the most boring one is probably the most accurate.

The people who worry that the iPad will bring about a dystopian future for home computing keep forgetting something: for the rest of humanity, their ideal world of perfectly hackable machines is already a dystopian nightmare. It’s a world in which nothing works without spending hours setting it up, in which basic features are missing while the manual lists thousands of irrelevant options, in which a million hardware extensions are available for their machine, but none of them help to solve a single one of their day-to-day problems. While being something of a hacker myself, I feel that the hacker’s vision of totally open computing probably should become a niche market, in much the same way that chemistry sets represent a niche market. The fact that not every person has a set of tools in his house that, by default, allows him to conduct arbitrary chemistry experiments has not substantially slowed down the progress of chemistry from what I can tell. The arrival of a world in which the most popular computers are closed to arbitrary hardware extensions and all applications are required to run within a sandbox probably won’t slow down the progress of personal computing much either.

Hackers of the world, your priorities are not simply different from the average user’s: they often represent a direct attack on the average user’s preferences. You keep asserting that you have the normal person’s interests in mind, but I think you’re often simply concealing your own self-interest underneath politicized rhetoric about freedom and openness.

Relevant to my earlier post about the relationship between direct democracy and laws prohibiting gay marriage, Pew Research just published poll data showing that a majority of Americans disapprove of same-sex marriage.

H. G. Wells famously said that, “statistical thinking will one day be as necessary for efficient citizenship as the ability to read and write.” I think we’re getting closer to that day: even the Supreme Court of the United States plans to start using the word ‘orthogonal’ colloquially.

Given the recent votes on same-sex marriage in New Jersey and Portugal, I wanted to test a seemingly innocuous claim that touches upon very broad issues in political theory: does the degree of directness of a “democratic” vote predict whether the vote will promote or prohibit same-sex marriage? Naively, it seemed clear to me that this was true: every single direct vote in my memory has outlawed same-sex marriage, while the only decisions that have allowed same-sex marriage have originated among elected representatives or unelected judges.

I decided to compile some rough data to test this idea using Wikipedia’s articles on the topic. I took what seemed to be the non-redundant decisions since the 1990’s from the following three articles:

1. Same-Sex Marriage in the Netherlands
2. Same-Sex Marriage in Spain
3. Same-Sex Marriage in the United States

And I arrived at this table:

You can find a CSV file with this data at the GitHub repository where I’ve stored all of the analyses I’ve completed so far. I’d love for people to add more data or contribute some visualizations of the patterns I’m pointing out. And, as always, I’d love to hear criticisms or suggestions about my approach.

The summary statistics provide pretty clear support for my intuition, though they’re rarely statistically significant, given the small sample sizes:

• 100% of court rulings promoted same-sex marriage.
• 80% of congressional decisions promoted same-sex marriage.
• 0% of direct votes promoted same-sex marriage.

For me the conclusion to be drawn from all of this is not that our respect for the will of the majority should compel us to prohibit same-sex marriage, but rather that we should be more hostile to direct democracy, because it is a powerful force for promoting intolerance in American society. I imagine this will sound un-American to many readers, but I think it’s quite clear that this sentiment was foundational for the American experiment in government. You see this distrust of direct democracy repeatedly reiterated in the Federalist Papers, especially in Federalist Number 10. Indeed, one can argue that the major purpose of our Constitution is to limit the ability of direct democracy to harm minorities within our society.

In further support of the claim that a lot of deaths are partly self-induced, here’s a fascinating piece by Wired on the extraordinary rise in the percent of deaths among the young caused by their own poor decisions.

It’s remarkable that, for the young, modern science has already made the world so safe that humanity, rather than nature, is now responsible for a majority of its own suffering.

My apologies for the long break between posts. Before the end of this week I’ll return to my series of posts on image processing in R. In the intervening time, I’ve finished a piece of code that I’d like to officially release to the public.

The code in question is a very minimal object-relational mapper written entirely in R. As of today, you can find the code on GitHub here. If you’re not familiar with using an ORM, I’d suggest reading a bit about how to use ActiveRecord, the ORM I’ve tried to emulate.

As it stands, the code I have is able to connect to a MySQL database and extract information about a specified table using MySQL’s SHOW COLUMNS. The code then automatically builds up R code that will map the rows of the table onto R objects. This code can subsequently be eval’d or cat’d to a file.

Let me note right off the bat that the code I’ve written is fairly heinous stylistically: because my understanding of R metaprogramming is rather limited, I’ve produced what amounts to the ugliest sort of code generation tool. If you check out the source in orm.R, you’ll quickly see what I mean: the output R code is generated as a string using a long series of paste() operations, many of which embed paste() operations inside of themselves. In the future, I would like to rewrite the code in a clearer fashion, possibly constructing expressions and parse trees directly, rather than using an intermediate layer of code as string.

That said, I want to release the code so that I can get feedback on my approach. To help you give me feedback, here’s a quick walkthrough of how’d you use the current code:

(1) Download the code from its GitHub home. Place the files you get into a directory of your choice. Let’s assume you use ~/r-ORM.

(2) Install the R YAML and MySQL packages if you don’t already have them on your system:

1
2

install.packages('yaml')
install.packages('RMySQL')

(3) Create a test database in MySQL called `sample_database`:

1

CREATE DATABASE `sample_database`;

(4) Give permissions on this test database to `sample_user`:

1
2
3

GRANT ALL ON `sample_database`.* TO 'sample_user'@'localhost' IDENTIFIED BY 'sample_password';
 
FLUSH PRIVILEGES;

(5) Create a test table called `users` in `sample_database`:

1
2
3
4
5
6
7
8

USE `sample_database`;
 
CREATE TABLE `users` (
  `user_id` INT(11) NOT NULL AUTO_INCREMENT,
  `name` VARCHAR(255) NOT NULL,
  `password` VARCHAR(255) NOT NULL,
  PRIMARY KEY (`user_id`)
) ENGINE=MyISAM DEFAULT CHARSET=latin1;

(6) Edit the database.yml file if you didn’t use the database name, user name or password I just suggested. You’ll also need to edit it if you’re not working on localhost.

(7) Open up an R interpreter. Set your working directory and source orm.R:

1
2

setwd('~/r-ORM')
source('orm.R')

(8) Build code for our user objects using orm.build.model('user'):

1
2

code <- orm.build.model('user')
cat(code)

At this point, you can review the code that’s been generated to see whether it should work on your system. If it will, you can eval it now:

1

eval(parse(text = code))

With that done, you should have a working set of functions that handle creating, finding, manipulating and deleting R objects that are serialized to the database. To test out the resulting model, let’s start by creating a user object. In general, an object of class foo will be created using an auto-generated function called create.foo:

1

user <- create.user()

Calling this functions builds a user object in memory. Nothing is in the database so far. To see the object, type user at the command line:

1

user

Now you can edit this user to make it a real piece of data. The columns of the database table are mapped onto object attributes with appropriate getter and setter methods, like so:

1
2

name(user) <- 'test_user'
password(user) <- 'test_password'

Once the user object is worth keeping, we store it in the database using store:

1

user <- store(user)

This stores the user object in the database. To get the ID’s edited correctly, you have to perform the assignment as indicated above. In the future I may change this.

After storing something, you might want to retrieve it later. Since we know that we just created the first user object, we can get it again by using a find.user call:

1

user <- find.user(1)

In general, you can find objects of class foo by calling find.foo(). If you provide an integer as an input, you’ll get the object with that ID. If you provide the string 'all', you’ll get a list containing all of the objects in your database. Other inputs produce an error.

You can see that we got the correct object using the getter methods:

1
2
3

user.id(user) == 1
name(user) == 'test_user'
password(user) == 'test_password'

We can edit it again to see that updating rows of the database works as expected:

1
2

name(user) <- 'new_user'
store(user)

Finally, now that we’re done with it, we can delete it from the database:

1

delete(user)

That, in a nutshell, is the use of this ORM solution. If you have any questions, please let me know: I’ll be happy to answer them.

I should note that there are many weaknesses of functionality with the current code. Notably:

1. If you want a ‘foo’ object, it must be in a ‘foos’ table. This is inspired by ActiveRecord’s default pluralization rule.
2. There’s no type checking. You can try to set string values on numeric columns; you’ll just get strange errors as a result.
3. There’s no checking whether you’re inserting a NULL into a column that won’t allow NULL’s.
4. The system might clobber an existing generic method or function if there’s a naming conflict. I’d recommend writing everything to disk before using this code until you’re familiar with the results.
5. The default objects are completely blank: the database default values are not used to initialize any attibutes of a new object.
6. The SQL generated is not sanitized. You should not use this system to process end-user input.

That said, I find the system usable for my current needs. I hope it will help you, and I really hope that you’ll consider making suggestions for hoping to improve the code. Patches would be especially welcome.

I would also love feedback about the interface that’s shown to the end user of this code. Here are some outstanding questions I’d like to hear from people about:

1. Are function names like orm.build.model() memorable?
2. Does the syntax of calls on the resulting objects seem sufficiently R-like?
3. Would it be useful to break out my minimal database abstraction layer into a separate package? Is there already a package to do this that I should be using instead?
4. How can I avoid the “code as string” approach I’m taking?
5. Do you have ideas about the ideal user syntax for implementing relations across tables, i.e. has_one or has_many relationships?
6. Should a caching mechanism be designed?
7. Should the database queries be minimized using some sort of pooling?
8. Should transactions be used?

This post is a continuation of my series dealing with matrix operations for image processing. My next goal is to demonstrate the construction of simple low-pass and high-pass spatial frequency filters in R. It’s easy enough to construct simple versions of these filters in R using the Fast Fourier Transform (also known as the FFT), but, because the FFT is a slightly complicated tool, I’m going to build up to using it progressively over a few posts.

For starters, I want to review the use of complex numbers in R. As always, if you’re interested in reviewing the mathematics of complex numbers, I’d start by browsing references online. The tutorial here seemed good to me at first glance, though I can’t claim to have read it through.

Because complex numbers are implemented in the “base” package, it’s very easy to start working with them. To construct the complex number x + iy, you use complex:

1
2
3
4
5
6
7

x <- 1
y <- 1
 
z <- complex(real = x, imaginary = y)
 
z
# [1] 1+1i

It’s conventional in mathematics to use z to refer to a complex number, so I’ll continue on with that tradition.

As always occurs with mathematical data types in R, you can convert other objects to class “complex” using as.complex:

1
2

as.complex(-1)
# [1] -1+0i

And you can test that an object is complex using is.complex:

1
2

is.complex(as.complex(-1))
# [1] TRUE

Beyond those standard operations, there are five essential mathematical operations you’d want to use on complex numbers.

First off, you want to be able to extract the real and imaginary components of a complex number. You can do this using Re and Im respectively:

1
2
3
4
5
6
7

z <- complex(real = 2, imaginary = 1)
 
Re(z)
# [1] 2
 
Im(z)
# [1] 1

The (x, y) representation of numbers is easier to understand at first, but a polar coordinates representation is often more practical. You can get the relevant components of this representation by finding the modulus and complex argument of a complex number. In R, you would use Mod and Arg:

1
2
3
4
5
6
7
8
9
10

z <- complex(real = 0, imaginary = 1)
 
Mod(z)
# [1] 1
 
Arg(z)
# [1] 1.570796
 
pi / 2
# [1] 1.570796

This corresponds to the intuition that i should be at a distance 1 from the origin and an angle of pi / 2.

Finally, you’ll want to be able to take the complex conjugate of a complex number; to do that in R, you can use Conj:

1
2
3
4
5

Conj(z)
# [1] 0-1i
 
Mod(z) == z * Conj(z)
# [1] TRUE

As you can see, the modulus of z equals z times the conjugate of z, which is exactly what you expect.

Now, historically, complex numbers were invented so that you could find the square root of negative numbers. By default, sqrt does not return a complex number when you ask for the square root of a negative number. Instead, it produces a NaN error, as you can see below:

1
2
3
4

sqrt(-1)
# [1] NaN
# Warning message:
# In sqrt(-1) : NaNs produced

To get the complex square root, you need to cast your negative number as a complex number using as.complex before applying sqrt:

1
2
3

sqrt(as.complex(-1))
 
# [1] 0+1i

In practice, complex numbers can be used to solve a huge range of problems, not least of which is the efficient representation of the FFT of an input vector. For fun, I thought I’d show a simple application: building a fractal using the Mandelbrot set.

To quote Wikipedia,

Mathematically the Mandelbrot set can be defined as the set of complex values of c for which the orbit of 0 under iteration of the complex quadratic polynomial z_n+1 = z_n² + c remains bounded. That is, a complex number, c, is in the Mandelbrot set if, when starting with z₀ = 0 and applying the iteration repeatedly, the absolute value of z_n never exceeds a certain number (that number depends on c) however large n gets.

We can translate this definition into R pretty easily by making certain assumptions about the exactness we want in our results. For my purposes, I’ll say that a number is in the Mandelbrot set if, after 100 iterations of the Mandelbrot algorithm, it’s never once had a modulus greater than 1,000,000. As you’ll see from the image I produced, this assumption works out pretty well, though it takes some real number crunching to get results out of it for reasonable sized data sets.

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

in.mandelbrot.set <- function(c, iterations = 100, bound = 1000000)
{
  z <- 0
 
  for (i in 1:iterations)
  {
    z <- z ** 2 + c
    if (Mod(z) > bound)
    {
      return(FALSE)
    }
  }
 
  return(TRUE)
}

Once we have this algorithm, we can easily generate an image of the Mandelbrot set by producing a matrix of complex numbers and depicting their inclusion in the set using image:

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

resolution <- 0.001
 
sequence <- seq(-1, 1, by = resolution)
 
m <- matrix(nrow = length(sequence), ncol = length(sequence))
 
for (x in sequence)
{
  for (y in sequence)
  {
    mandelbrot <- in.mandelbrot.set(complex(real = x, imaginary = y))
    m[(x + resolution + 1) / resolution, (y + resolution + 1) / resolution] <- mandelbrot
  }
}
 
png('mandelbrot.png')
image(m)
dev.off()

And here’s the resulting image:

I’m not quite sure what the source of the lines in my image is. I’ll have to report back on that.

I’ve said it before: hate speech laws, whenever they exist, will always be exploited to suppress speech that’s politically charged. Italy’s proving my point this week, as it contemplates Web restrictions after the attack on Berlusconi.

Over the next few posts, I’m going to be reviewing the use of R to implement the most commonly used matrix techniques for image manipulation. The code will be surprisingly simple to understand, because the real magic behind these techniques lies in the mathematics that R provides an abstract interface to. To start, I’m going to review the singular value decomposition of a matrix, also known as the SVD.

If you’d like to understand the magic behind the scenes here, you can read about the mathematical definition of the SVD on Wikipedia or you can watch Gilbert Strang’s lectures. If your command of linear algebra isn’t great, I suspect that the mathematical details of the SVD will be totally opaque to you.

Fortunately, the details of the underlying mathematics aren’t at all necessary for appreciating the SVD’s usefulness for manipulating images: it is arguably the best way to compress matrices into smaller representations that contain the essential information from the original matrices. The easiest way to understand this is to see it in action, so I’m going to show how the SVD allows for any degree of compression of an image represented as a real-valued matrix.

Obviously, the first thing we have to do is to represent our example image as a matrix with real-valued entries. To do this, we’ll need to use a simple technique for loading an image into R that I found on the Quantitative Ecology blog. The first thing to do is use ImageMagick to transform a TIFF file into a PPM file:

1
2
3

#!/bin/bash
 
convert face.tiff image.ppm

Then we can use the “pixmap” package from CRAN to read the PPM image into R as an image object containing three matrices: one for the red components, one for the green components, and one for the blue components. The exact code I used was this:

1
2
3
4
5
6
7

library('pixmap')
 
image <- read.pnm('image.ppm')
 
red.matrix <- matrix(image@red, nrow = image@size[1], ncol = image@size[2])
green.matrix <- matrix(image@green, nrow = image@size[1], ncol = image@size[2])
blue.matrix <- matrix(image@blue, nrow = image@size[1], ncol = image@size[2])

Before we go any further, we should look at the image we’ll be playing with. Rather than work with the composite of the three color matrices, I found it easier to work with only a single color. For this specific image, I found that the green matrix looked best, so I’ll use only that. At the risk of seeming self-indulgent, I’m using my own Gravatar icon as an example image because faces are particularly easy to recognize under the strange color scheme I’ll be using, which I chose because it really simplifies the visualization code.

To visualize the matrix as an image, I discovered two very helpful graphical functions, image and heat.colors. image depicts the value of a matrix using a color scheme you specify; heat.colors produces a set of colors that could be used in a heat map with a granularity you specify as the first argument. Using these functions, we can see our raw input matrix:

1

image(green.matrix, col = heat.colors(255))

The images generated by image are at an unfortunate angle for viewing, so I edited them post hoc with a single 90º clockwise rotation in Preview to give this sort of image:

If you want to play with the original PPM image file, you can get it here.

Clearly, we’ve got a recognizable face here. Feeling good about that, we can start to play with the SVD of this matrix. The R function to generate the SVD of a matrix is simply svd:

1

green.matrix.svd <- svd(green.matrix)

svd returns a list with three entries: d, u and v. d is a vector of the singular values ordered by decreasing size; u and v are matrices that are combined with a diagonal matrix form of d to give the original input matrix. To simplify things, I’ll set up variables to contain the value of these elements of the SVD output list:

1
2
3

d <- green.matrix.svd$d
u <- green.matrix.svd$u
v <- green.matrix.svd$v

We can check that the SVD works by multiplying the relevant matrices:

1
2
3
4

green.matrix.reconstruction <- u %*% diag(d) %*% t(v)
 
mean(green.matrix - green.matrix.reconstruction)
# [1] -5.184204e-16

If this is the SVD by itself, it’s not clear that it’s useful computationally: we’ve just rewritten one matrix as the product of three matrices and accumulated a small bit of error along the way. The real utility of the SVD lies in the singular values: they represent, in decreasing order, the most important information about the original matrix.

To see this, you can shrink the input matrices and produce a compressed form of the matrix. For instance, my original matrix is a 212 by 201 image. Using the SVD, we can represent it as the product of a 212 by 2 matrix, a 2 by 2 matrix, and a 2 by 201 matrix. This massively reduces the amount of information we’re storing, but you can already see the outline of our input image in the results:

1
2
3

green.matrix.compressed <- u[,1:i] %*% diag(d[1:i]) %*% t(v[,1:i])
 
image(green.matrix.compressed, col = heat.colors(255))

To get a feel for how this works, we can iterate over the number of singular values we include in our compressed form:

1
2
3
4
5
6
7

for (i in c(3, 4, 5, 10, 20, 30))
{
	green.matrix.compressed <- u[,1:i] %*% diag(d[1:i]) %*% t(v[,1:i])
	png(paste('images/Green Image SVD Compressed ', i, '.png', sep = ''))
	image(green.matrix.compressed, col = heat.colors(255))
	dev.off()
}

This produces the following images:

3 Singular Values

4 Singular Values

5 Singular Values

10 Singular Values

20 Singular Values

30 Singular Values

As you can see, the images get better quite quickly. Indeed, this last image looks great, despite the large reduction in the number of entries we’re keeping from the original matrix. I’d say that’s a pretty compelling case for the value of the SVD. It would be easy to combine the compressed matrices for red, green and blue to get a compressed version of our input image using its original color scheme. I’ll leave that as an exercise for anyone interested in playing with this technique on their own.

Of course, the SVD has tons of other uses, but this simple hack for image compression struck me as pretty interesting, as well as being remarkably simple to implement in R.