John Myles White: Die Sudelbücher

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.

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.

Lately, I’ve been running a series of fMRI experiments on visual perception. In the interests of understanding the underlying properties of the images I’m using as stimuli, I’ve been trying to learn more about the matrix transformations commonly used for image compression and image manipulation. Thankfully, R provides simple-to-use implementations for all of the matrix operations I wanted to play around with, so it’s been quite easy to get started. For the next few posts, I thought that I’d review the standard matrix techniques for image compression and editing, giving full examples of their implementation in R and demonstrations of their real-world value.

Before I start, I should make sure that you’re familiar with basic matrix operations in R. I imagine almost every R user knows a little bit about matrix algebra and probably knows the basics of using R to perform matrix algebra, but here’s a quick review to make sure I don’t leave anyone in the dark:

Building Matrices

You can build a matrix in R using the matrix function:

1
2
3
4
5
6

m <- matrix(c(1, 0, 0, 1), nrow = 2, ncol = 2, byrow = TRUE)
 
m
#      [,1] [,2]
# [1,]    1    0
# [2,]    0    1

You can test whether an item you’ve been given is a matrix using is.matrix and you can convert appropriate objects to matrices using as.matrix:

1
2
3
4
5
6
7
8
9

m <- matrix(c(1, 0, 0, 1), nrow = 2, ncol = 2, byrow = TRUE)
 
is.matrix(m)
# TRUE
 
as.matrix(data.frame(x = c(1, 0), y = c(0, 1)))
#      x y
# [1,] 1 0
# [2,] 0 1

Diagonal Matrices

The matrix I’ve been building in the examples above is a diagonal matrix, so you could also construct it using diag:

1
2
3
4
5
6

m <- diag(nrow = 2, ncol = 2)
 
m
#      [,1] [,2]
# [1,]    1    0
# [2,]    0    1

In general, you can generate the n by n identity matrix as:

1
2
3

n <- 10
 
diag(nrow = n, ncol = n)

I’ll be exploiting a more interesting use of diag in my next post, so let’s see how you can build matrices other than the identity with diag by specifying a vector of entries along the diagonal:

1
2
3
4
5
6

m <- diag(c(2, 1), nrow = 2, ncol = 2)
 
m
#      [,1] [,2]
# [1,]    2    0
# [2,]    0    1

This form of diag turns out to be extremely useful, as you’ll see once I cover the SVD’s syntax in R.

Matrix Algebra: Addition, Scalar Multiplication, Matrix Multiplication

The three core operations that can be performed on matrices are addition, scalar multiplication and matrix multiplication. These are easy to work with in R:

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

m <- matrix(c(0, 2, 1, 0), nrow = 2, ncol = 2, byrow = TRUE)
 
m + m
#      [,1] [,2]
# [1,]    0    4
# [2,]    2    0
 
2 * m
#      [,1] [,2]
# [1,]    0    4
# [2,]    2    0
 
m %*% m
#      [,1] [,2]
# [1,]    2    0
# [2,]    0    2

Be careful with the * operator: it does not perform matrix multiplication, but rather an entry-wise multiplication:

1
2
3
4
5

m * m
 
#      [,1] [,2]
# [1,]    0    4
# [2,]    1    0

Matrix Transposes and Inverses

The next few matrix operations are a little more complex: transposition and inversion. Transposition is the easier of the two. To get the transpose of a matrix, you simply call the t function:

1
2
3
4

t(m)
#      [,1] [,2]
# [1,]    0    1
# [2,]    2    0

In contrast, inversion is a little more complex, partly because the function you’d want to use has a non-obvious name: solve.

1
2
3
4
5

solve(m)
 
#      [,1] [,2]
# [1,]  0.0    1
# [2,]  0.5    0

The reason that solve is called solve is that it’s a general purpose function you can use to solve matrix equations without wasting time computing the full inverse, which is often inefficient. If you want to know more about the computational efficiency issues, you should look into the ideas behind the even faster variant, qr.solve.

Now, you probably know this already, but the definition of a matrix’s inverse is that the product of the matrix and its inverse is the identity matrix, if the inverse exists. I always find this a good way to make sure that I’m correctly computing the inverse of a matrix:

1
2
3
4
5

solve(m) %*% m == diag(nrow = nrow(m), ncol = ncol(m))
 
#      [,1] [,2]
# [1,] TRUE TRUE
# [2,] TRUE TRUE

Further Matrix Operations

The car package defines an inv function, which is simply a new name for solve:

1
2
3
4
5
6

library('car')
inv
 
# function (x) 
# solve(x)
# <environment: namespace:car>

I think this is pretty clever, though I’m loathe to import a package just for this one snippet.

More interestingly, the MASS package defines a ginv function, which gives the matrix pseudoinverse, a generalization of matrix inversion that works for all matrices:

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

library('MASS')
 
m <- matrix(c(1, 1, 2, 2), nrow = 2, ncol = 2)
 
solve(m)
# Error in solve.default(m) : 
#  Lapack routine dgesv: system is exactly singular
 
ginv(m)
#      [,1] [,2]
# [1,]  0.1  0.1
# [2,]  0.2  0.2

In practice, you can usually get around using the pseudoinverse, but it’s nice to know that it’s at hand all the time.

Eigenvalues and Eigenvectors

Eigenvectors are surely the bane of every starting student of linear algebra, though their considerable power to simplify problems makes them the darling of every applied mathematician. Thankfully, R makes it easy to get these for every matrix:

1
2
3
4
5
6
7
8
9
10

m <- diag(nrow = 2, ncol = 2)
 
eigen(m)
# $values
# [1] 1 1
#
# $vectors
#      [,1] [,2]
# [1,]    0   -1
# [2,]    1    0

Matrix Metadata

Last, but not least, you can get metadata about the shape of a matrix using dim, nrow and ncol:

1
2
3
4
5
6
7
8
9
10

m <- diag(nrow = 2, ncol = 2)
 
dim(m)
# [1] 2 2
 
nrow(m)
# [1] 2
 
ncol(m)
# [1] 2

Hopefully those operations were already familiar to any readers out there, as I doubt that they’ll be clear from this short explanation without prior knowledge. I just felt compelled to review them before using any of them in my next set of posts. Tomorrow, I’ll start going through more interesting matrix algorithms in R, beginning with the SVD.

With a little guidance from the indefatigable Hadley Wickham, I figured out today how to implement the setter methods that were missing from my example user class. To review, let’s rebuild the getter methods for my user object:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

user <- list(id = 1,
             password = '286755fad04869ca523320acce0dc6a4',
             email = 'jmw@johnmyleswhite.com')
 
class(user) <- 'user'
 
id <- function(x)
{
  UseMethod('id', x)
}
 
id.user <- function(user)
{
  return(user[['id']])
}
 
password <- function(x)
{
  UseMethod('password', x)
}
 
password.user <- function(user)
{
  return(user[['password']])
}
 
email <- function(x)
{
  UseMethod('email', x)
}
 
email.user <- function(user)
{
  return(user[['email']])
}

With these working, my goal was to write setter methods that would work like the ideal code below:

1
2
3
4
5
6

id(user) <- 2
 
if (id(user) == 2)
{
  print("Succeeded in editing the user's id attribute.")
}

Defining this sort of setter method turned out to require a little effort. I wasted quite a lot of time walking down dead ends, but, thankfully, Hadley gave me the exact piece of information I was missing early this morning. The dead ends were interesting in themselves, though, so I’ll review them in another post tomorrow.

For now I’ll just explain the correct implementation for my desired setter methods, which you can see below:

1
2
3
4
5
6
7
8
9
10

'id<-' <- function(x, value)
{
  UseMethod('id<-', x)
}
 
'id<-.user' <- function(x, value)
{
  x[['id']] <- value
  return(x)
}

There are three important ideas here:

1. The generic setter method for id should be called 'id<-', which must always be quoted to prevent parsing errors.
2. The class level setter method should be called 'id<-.user', which isn't surprising, though I had imagined it might be called 'id.user<-' at one point.
3. The setter method should make a change to a copy of the original object and return the edited copy to the caller. R handles the assignment of this edited return value to the original object behind the scenes. In fact, not using this return value approach will make otherwise plausible looking code fail to edit your object correctly. I stumbled on that for quite a while before I was shown the light.

With this in mind, a final user class looks like this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

user <- list(id = 1,
             password = '286755fad04869ca523320acce0dc6a4',
             email = 'jmw@johnmyleswhite.com')
 
class(user) <- 'user'
 
id <- function(x)
{
  UseMethod('id', x)
}
 
id.user <- function(user)
{
  return(user[['id']])
}
 
password <- function(x)
{
  UseMethod('password', x)
}
 
password.user <- function(user)
{
  return(user[['password']])
}
 
email <- function(x)
{
  UseMethod('email', x)
}
 
email.user <- function(user)
{
  return(user[['email']])
}
 
'id<-' <- function(x, value)
{
  UseMethod('id<-', x)
}
 
'id<-.user' <- function(x, value)
{
  x[['id']] <- value
  return(x)
}
 
'password<-' <- function(x, value)
{
  UseMethod('password<-', x)
}
 
'password<-.user' <- function(x, value)
{
  x[['password']] <- value
  return(x)
}
 
'email<-' <- function(x, value)
{
  UseMethod('email<-', x)
}
 
'email<-.user' <- function(x, value)
{
  x[['email']] <- value
  return(x)
}

Until recently, I’ve never had any reason to learn how to define my own classes in R. Having learned this week, I was surprised to find out how easy it is to start implementing classes in R. If you know nothing about creating classes and class methods in R, here’s a very quick overview of the three core ideas behind R’s object system. If you already know the basics of defining classes in R, I’d suggest skipping this post, as it’s not likely to be of much value to you.

The first thing you should know is that the S3 object-oriented system feels like a hack that was put on top of the original S language. If you’re familiar with Perl, you’ll feel right at home with the basic ideas behind this approach; you’ll probably also appreciate the conceptual simplicity of the results.

With that in mind, you should be aware that any normal piece of R code is filled with objects, because the base packages and most well-designed user packages use lots of custom classes. You can see this for yourself by starting to use class to perform introspection on some of the items in your programs:

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

class(1)
# [1] "numeric"
 
class('a')
# [1] "character"
 
x <- 1:10
y <- 5 * x + rnorm(length(x), 0, 1)
 
fit <- lm(y ~ x)
 
class(fit)
# [1] "lm"

More interesting than looking at the classes of existing objects is defining new classes of your own. Thankfully, class, like names and many other functions in R, can be assigned to directly, like so:

1
2
3
4
5
6

a <- 1
 
class(a) <- 'octonion'
 
class(a)
# [1] "octonion"

If you’re familiar with Perl, you can think of this assignment to class as analogous to using bless. If you’re not a Perl user, hopefully this approach to defining classes still makes sense to you, because it’s so simple: the only thing that decides the class of an object is the value you’ve set for the class attribute. To convince yourself that the class of an object is just an attribute, you can use attr:

1

attr(a, 'class')

This simple metadata driven approach to building objects is the core of the S3 object system, as far as I can tell.

Of course, setting class by itself is pretty useless: you want to be able to define class methods. That’s where the other two main ideas come in. The first trick is to use generic methods to get polymorphism out of the language; the second trick is to define methods on your classes using a simple naming convention. To see how this works, let’s use an example that should be familiar to anyone who’s ever built a database-backed website.

Instead of my would-be octonion class, we’ll consider a user object. To define an object of class ‘user’, we can do the following:

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

user <- list(id = 1,
             password = '41bfe136a536b7749104415bc364df2e',
             email = 'jmw@johnmyleswhite.com')
 
class(user) <- 'user'
 
user
 
#$id
#[1] 1
#
#$password
#[1] "41bfe136a536b7749104415bc364df2e"
#
#$email
#[1] "jmw@johnmyleswhite.com"
#
#attr(,"class")
#[1] "user"

The first thing that comes to mind after building this object is that we should define getter and setter methods for accessing and modifying the contents of this user object. For example, we want a way to get the id attribute for our object. We’d also like our approach to generalize to other objects with an id attribute. Ideally, we should be able to write something like this:

1

id(user)

Now, it is obviously possible to define an id function that operates only on user objects, but that’s not the right approach if we also want to have another sort of object that would have an id method as well. This polymorphism concern is the problem that generic methods and specialized naming conventions solve.

First, we are going to define a specialized id method that only operates on user objects. Because it will only work on ‘user’ objects, we’ll call it id.user:

1
2
3
4
5
6
7

id.user <- function(user.object)
{
	return(user.object[['id']])
}
 
id.user(user)
# [1] 1

Given this, we can then define a generic method that will operate on objects of many classes and reroute our general function calls to the correct class-level method:

1
2
3
4
5
6
7

id <- function(object)
{
	UseMethod('id', object)
}
 
id(user)
# [1] 1

Here UseMethod just searches for an id.user function, finds it and then calls it with user as its argument. If we had a ‘profile’ object, id would search for an id.profile function and call that. You can see this by trying id on a variable whose class we set to be profile using class.

1
2
3
4
5
6
7
8

profile <- list(id = 1)
 
class(profile) <- 'profile'
 
id(profile)
 
# Error in UseMethod("id", user) : 
#  no applicable method for 'id' applied to an object of class "profile"

With these ideas in mind, it’s easy to do something similar for the rest of our attributes as well:

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

password.user <- function(user.object)
{
	return(user.object[['password']])
}
 
password <- function(user)
{
	UseMethod('password', user)
}
 
email.user <- function(user.object)
{
	return(user.object[['email']])
}
 
email <- function(user)
{
	UseMethod('email', user)
}

Defining setter methods is a little more tricky. We could use a trick like the eval hacks I used to implement push and pop recently, but Hadley Wickham wisely chastened me for doing that. I’m still trying to decide what’s the best approach.

As I see it, there are least three possible method call styles you might define:

1
2
3
4
5

user <- password(user, 'new_value')
 
password(user, 'new_value')
 
password(user) <- 'new_value'

I would really like to start implementing this last sort, but I haven’t figured out how to yet. The others can be implemented using the generic methods approach I just outlined above, along with some ugly calls to eval for the second call style. To implement the first, simply do something like this:

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

password.user <- function(user, new.value = NULL)
{
  if (is.null(new.value))
  {
    return(user[['password']])
  }
  else
  {
    user[['password']] <- new.value
    return(user)
  }
}
 
password <- function(user, new_value = NULL)
{
  UseMethod('password', user, new_value)
}

If you have suggestions on how to implement the third approach, please let me know. Also, I am not at all happy with the use of NULL in my current setter implementation, because that makes it impossible to set the value of an attribute to NULL. If people have suggestions, I’d very much appreciate them. Does R implement arity-specific function definitions?