What's Watson Developer Cloud

A few times in the last year or so an IBM sales-person has come around to the company I work at spruiking products branded with the IBM 'Watson' label. To be frank, I came out of these meetings without a clear idea of what, exactly, the Watson suite of products was, or how it could be used in our business. The problem was that there were too many products sitting under the Watson brand, and the sales rep kept changing between different things.

So, I've just spent a couple of days in Melbourne, getting a better handle on what is available, what it can be used for, what it costs, etc. And here is a quick summary:

What is Watson?

I still can't definitively answer this. It seems to just be a brand/label for a range of products/APIs, focussed an analytics of some sort (speech analysis, text analysis, image analysis, sentiment analysis, etc). There are a few things on offer from IBM, and I think it helps to separate those:

1. Services: IBM has a range of analytical services (language processing, image recognition, and so on), all exposed through APIs, and you pay per API request to use these. So you can submit (say) an image to the Watson image classification service and it tells you the probability that the image contains a cat (say). More on this in a moment.
2. Application hosting: you can build an application that uses the IBM services, and IBM can also host that application for you. I'm not going to say more about this, as it is of less interest to me right now.

I have to say that the communication and documentation around all this from IBM is not good, and to start with it's all a bit bewildering as you work out how tightly coupled the app hosting is to the services (do you want an app, a virtual server? a container? a service? do you want to 'bind' your hosted app to your service? Too too many things and no coherent presentation of what they are and how they work together.... for example, it's still not clear to me when you would bind a service to an app and when you wouldn't.... it seems that you can use a service from an un-bound application, so maybe it's just a performance thing?). But anyway, I digress... from now on I'm just going to talk about the services, not the apps.

What Services are available? What do they do?

Below, I list the services that are available under the Watson brand, as at time of writing. There is an API that exposes each of the services, and the model seems to be that you are changed per API transaction, for each service. So for instance, there is an image classification service that you can train with some data, and then submit images to that service to get classification results. Or you can submit free text to the tone analysis service and it will return you information about the 'tone' of the language (aggressive, friendly, etc). Each such submitted sample would cost you an API call, which costs between $0.2 and $0.05 (depending on the volume of calls, you get a discount for volume). Pricing varies by service. You can get up to date info on pricing if you dig down from here.

Here are the services on offer, each of which you can call via a simple API (more on this in a second). The description next to each is the IBM-provided one. I'll say a bit more about a few of them underneath.

Concept Expansion: Maps euphemisms or colloquial terms to more commonly understood phrases

Concept Insights: Explore the concepts behind your input, identifying associations beyond traditional text matching.

Dialog: Enable your application to use natural language to converse with users

Document Conversion: Converts a  HTML, PDF, or Microsoft Word™ document into a normalized HTML, plain text, or a set of JSON-formatted Answer units.

Language Translation: Translate text from one language to another for specific domains.

Natural Language Classifier: Natural Language Classifier performs natural language classification on question texts. A user would be able to train their data and the predict the appropriate class for a input question.

Personality Insights: The Watson Personality Insights derives insights from transactional and social media data to identify psychological traits

Relationship Extraction: Intelligently finds relationships between sentences components (nouns, verbs, subjects, objects, etc.)

Retrieve and Rank: Add machine learning enhanced search capabilities to your application

Speech To Text: Low-latency, streaming transcription

Text to Speech: Synthesizes natural-sounding speech from text.

Tone Analyzer: It helps people detect, understand and revise the language tones of emotions, social propensities and writing styles from their writings.

Tradeoff Analytics: Helps make better choices under multiple conflicting goals. Combines smart visualization and recommendations for tradeoff exploration

Visual Recognition: Analyzes the visual content of images and videos to understand their content without requiring a textual description

Cognitive Commerce, Graph, and Insights: These are third-party provided. 

Some comments

Exposing analytics as a service could be great, but I don't think Watson quite gets you to where you want to be. Maybe it will soon, but not right now. For example, suppose you want to classify images into one of (say) 5 categories. Right now, you cannot do this with Watson except in a clumsy way -- the API appears to only support binary classification for images, so multi-class classification is not possible. Sure, you can build multiple binary classifiers but this is silly cludge to have to work through. There are some strange other restrictions too: you are limited to 15,000 training samples when you train the natural language classifier. Why? Training a deep network to do multi-class classification would certainly benefit from having more than 15,000 training samples, so why put this limit in there? Also, it's not clear to me what sort of training is going on under the hood. I know that putting things behind an API is meant to mask all that, but training neural network's can be quite fiddly and task specific... and I am really quite doubtful that the default network architecture, learning rate, batch size, etc are anywhere near optimal. It would be nice if IBM did hyper-parameter tuning for you behind the scenes too, but I'm pretty damn sure this is not happening -- you just get a trained classifier that may or may not have good settings for your particular problem.

In general, the Watson API offerings seem to occupy a space that will appeal to smaller organizations without a lot of analytics capacility, and may wow some managers who want 'insights' but don't care if those insights are accurate or not, or have a clear idea about what, exactly, they want to do with those insights. The prime example of this is the personality analysis API -- you can feed in free text (from politicians speeches, or a twitter handle, or wherever) and get back a 'personality model' of the person who wrote that text. The personality model is a bit like a horoscope though: almost anything it returns sounds half-plausible. So I can see managers saying 'great, this tells us which customers are open, or angry, or supportive, and this will let us target marketing to them!', but are the scored personality traits returned by the Watson API actually useful in targeted marketing? Who knows, but I'm certainly sceptical, and I'm scared that the richness of the outputs will be mis-taken by managers as something that is both accurate and useful, when it's far from clear that it is either.

   

Simplest possible Theano examples

This is another write-up of my steps learning Theano. My intention is to end up using an easier-to-use library such as keras, which sits on top of Theano, but before I do that I want to make sure I have a reasonable handle on Theano itself.

I started out with the logistic regression example available here, and then went on to the denoising auto-encoder example available here. However, being completely new to Theano, I thought the examples were not as simple as they could be, so I stripped the examples down to the bare bones (no classes, no functions, no bothering about proper weight initialization). I wanted to focus on understanding Theano itself, and the code in the existing example is a bit of a distraction from that. Others may find themselves in a similar position (wanted the simplest possible examples), and that's why I've written them up and put them here.

Example 1: Logistic Regression

There are a couple of logistic regression examples for Theano, but of those, I think this one is the simplest and best one to start on. It is a great place to start and easy to do in Theano while getting used to the main concepts. I've reproduced the code below, with some additional comments added by me.

import numpy
import theano
import theano.tensor as T
rng = numpy.random

# Create some dummy data. Here we create a dataset of size 400
# each example has 784 independent variables/columns, and a randomly
# assigned 'label' between 0 and 1 (inclusive), so we are doing binary
# logistic regression

N = 400   #400 samples
feats = 784 #each with 784 features
D = (rng.randn(N, feats), rng.randint(size=N, low=0, high=2))
training_steps = 10000

# Declare Theano symbolic variables
x = T.matrix("x")  # input matrix (of size (N, feats))
y = T.vector("y")  # output vector of length N
w = theano.shared(rng.randn(feats), name="w") # weights
b = theano.shared(0., name="b") # a single bias/constant term
print("Initial model:")
print(w.get_value())
print(b.get_value())

# Construct Theano expression graph
p_1 = 1 / (1 + T.exp(-T.dot(x, w) - b))   # Probability that target = 1
prediction = p_1 > 0.5                  # The prediction thresholded
xent = -y * T.log(p_1) - (1-y) * T.log(1-p_1) # Cross-entropy loss function
cost = xent.mean() + 0.01 * (w ** 2).sum() # The cost to minimize. Note this includes a L2 penalty on the weights

gw, gb = T.grad(cost, [w, b])  # Compute the gradient of cost wrt w/b 

# Compile. The 0.1 in here is the learning rate (basically how quickly we try to go downhill when doing optimisation
train = theano.function(
          inputs=[x,y],
          outputs=[prediction, xent],
          updates=((w, w - 0.1 * gw), (b, b - 0.1 * gb)))
predict = theano.function(inputs=[x], outputs=prediction)

# Train
for i in range(training_steps):
    pred, err = train(D[0], D[1])

print("Final model:")
print(w.get_value())
print(b.get_value())
print("target values for D:")
print(D[1])
print("prediction on D:")
print(predict(D[0]))

Example 2: Auto-encoder

I didn't like the structure of the auto-encoder tutorial available here, so here is the stripped down example I worked through to solidify my Theano understanding. The nice thing about the following code is that the whole thing works if you run each code segment in the order presented here.

Step 1: get MNIST data

We will train our auto-encoder on the MNIST data, so we are essentially going to create a single later network which creates some key features of the MNIST digits, and then used those features to reconstruct (as best as possible) the MNIST digits.

For the example below, I will assume that you have downloaded the MNIST data (mnist.pkl.gz) and saved it somewhere.

Step 2: import required libraries and write data import function

Lets first make all the imports that we are going to need, and then there is a function which will load in the MNIST data for us. This function does this for us.

print '... loading data'  
# Load the dataset  
f = gzip.open("/path/to/where/you/saved/mnist/mnist.pkl.gz", 'rb')  

train_set, valid_set, test_set = cPickle.load(f)  
f.close()  
data_x, data_y = train_set  
train_set_x = theano.shared(numpy.asarray(data_x,  
                  dtype=theano.config.floatX))   

Nothing exciting there -- we just have our training data set, which is a matrix with 50000 rows (i.e. images), and 784 columns, which are the values in a 28x28 image, flattened out into a single length 784 vector. The values are between 0 (black) and 1 (white).

Step 3: Define our auto-encoder 

Now we need to specify our auto-encoder structure. We have 784 inputs (28x28 images), and so must also have 784 outputs (because we are trying to reconstruct the source image). I've gone for 500 hidden units, for no particularly good reason.

Because we have 784 inputs/outputs and 500 hidden nodes, we need a weights matrix of size (784,500) between our inputs and hidden layer, and one of size (500,784) between our hidden and output layers. Below, these are called W1 and W2. We also have two bias vectors b1 and b2. b1 is a vector of biases for the hidden layer (so length 500), and b2 is a vector of biases for the output layer (so length 784).

You can get into a long debate about the best way to initialize these weights, but I'll just be lazy and initialize them to something random and reasonably close to 0. I'll return to a discussion of the 'correct' way to initialize weights later.

#attempt at a simple auto-encoder in Theano  
#  
#  
#We do this by mapping a 28x28 image (i.e. 784 inputs)   
#into a hidden layer, and then back out again to a visible (output)   
#layer  
n_visible=784  
n_hidden=500  
#These are our input images. Each row is of length 784  
#The number of rows will depend on the batch size of course  
x = T.matrix('x')  
#These are the weights between the input and the hidden layer. There must be   
#784*n_hidden of them  
W1 = theano.shared(value=(rng.rand(n_visible, n_hidden)-0.5)*0.1, name="W1")  
#This is the bias into the first hidden layer, so we need n_hidden of those  
b1 = theano.shared(value=rng.rand(n_hidden)*0.02-0.01, name="b1")  
#These are the weights between the hidden layer and the output layer  
W2 = theano.shared(value=(rng.rand(n_hidden, n_visible)-0.5)*0.1, name="W2")  
#These are the biases for the output layer  
b2 = theano.shared(value=rng.rand(n_visible)*0.02-0.01, name="b2")  
 

Now, in principle, W1 and W2 are completely separate matrices....but it's a common trick to tie the two together, so W1 = transpose(W2). The reasoning for this is intuitively appealing but not mathematically precise: W1 encodes the image, so it is nice to think of the transpose operation as 'unpacking'/decoding the image. This would work exactly if the hidden layer had a linear activation function and W1 was orthogonal, but this is not the case. Still, we cut the number of parameters we have to estimate in half if we tie the weights, so maybe it is not a bad idea to do that.

I pick a non-linear activation function for the hidden units: the rectified linear transformation. I pick a sigmoid function for the output node, so that the response is squashed back to between 0 and 1 (the same range as the input image).

#This is the output from the first layer 
L1output = T.nnet.relu(T.dot(x, W1) + b1)  
#This is our output from the second layer, if we use tied weights  
L2output = T.nnet.sigmoid(T.dot(L1output, W1.T) + b2)  
 

Now that we've specified the network, we need to say something about how to train the network. Just being lazy, I choose to minimize squared error, which is often an OK choice. Once I've made that choice, I just need to get the gradients of the network parameters (W1, b1, b2) w.r.t the loss, and that will be enough for me to train my network.

#OK, we've specified our weights/bias structure, now lets specify   
#how to do updates.   
#First, lets work out the loss at the output nodes, for input x  
#Here I just use squared loss  
loss = T.mean((L2output-x)**2)  
#compute gradients w.r.t input parameters  
W1grad, b1grad, b2grad = T.grad(loss, [W1, b1, b2])  

Now, we've specified the loss we want to minimise, and the gradients required to minimise that loss, but we need to put that together in a training function, which we specify like so:

#now specify our update/training step  
trainf = theano.function(  
    inputs=[x],  
    outputs=[L2output, loss],  
    updates=((W1,W1-0.1*W1grad),(b1,b1-0.1*b1grad),(b2,b2-0.1*b2grad))  
)  
predict = theano.function(inputs=[x], outputs=L2output)  

This just says that we update our weights and bias in the direction of the negative gradient. The 0.1 in the trainf function is just the learning rate -- higher values will mean faster learning (but maybe we get stuck in a local minima), lower values mean slower training but probably a better final model.

Now, that's pretty much it in terms of specifying the model and how to train it. Now we just need to do the actual training! This is pretty simple -- we break our data up into batched and do stochastic gradient decent on the batches:

#OK, that's it I believe! Now compile it and do the actual training  
#Before we do that, we just need to get Theano to 'evaluate' the data  
trainX = train_set_x.get_value()  
batch_size = 20  
for i in range(10):  
    batches = trainX.shape[0]/batch_size  
    errssofar = []  
    for batch in range(batches):      
        pred, err = trainf(trainX[batch*batch_size:(batch+1)*batch_size])  
        errssofar.append(err)  
        avgerr = numpy.mean(errssofar)  
        print("Iteration "+str(i)+" error on batch "+str(batch)+" is "+str(err)+". Average error across all batches this iteration is "+str(avgerr))  

Now that should get you through an initial 10 epochs of training. So now we'd like to see how well the auto-encoding has worked. Lets define some functions to let us look at the input images, the recovered output images, and the weights in the hidden neurons:

import Image  
#have a look at a particular image (input)  
def showInput(index):  
    Image.fromarray(numpy.asarray([int(item*255.99) for item in trainX[index]], dtype=numpy.uint8).reshape(28,28)).show()  
#have a look at a reconstructed output   
def showOutput(index):  
    Image.fromarray(numpy.asarray([int(item*255.99) for item in predict(trainX[index:index+1])[0]], dtype=numpy.uint8).reshape(28,28)).show()  
#have a look at some weights in the hidden layer  
def showHiddenUnit(index):  
    weights = [item for item in W1.get_value().T[index]]  
    minw, maxw = min(weights), max(weights)  
    weights = [int(255.99*(item-minw)/(maxw-minw)) for item in weights]  
    Image.fromarray(numpy.asarray(weights, dtype=numpy.uint8).reshape(28,28)).show()  

So lets now use these functions to look at one of our inputs, look at its reconstructed output, and then look at the weights:

showInput(0)  
showOutput(0)  
showHiddenUnit(0)  

Probably after only 10 iterations the reconstruction is not that good, so you might want to train a bit more before looking at the reconstruction again:

for i in range(10):  
    batches = trainX.shape[0]/batch_size  
    errssofar = []  
    for batch in range(batches):      
        pred, err = trainf(trainX[batch*batch_size:(batch+1)*batch_size])  
        errssofar.append(err)  
        avgerr = numpy.mean(errssofar)  
        print("Iteration "+str(i)+" error on batch "+str(batch)+" is "+str(err)+". Average error across all batches this iteration is "+str(avgerr)) 

Repeat this as many times as you like and see how the encoding gets a bit better with more training.

Initialisation of weights

Above, I just glossed over the initialisation of weights, but apparently this is quite an important thing, and there has been a lot of research into the effect of weights initialisation, and the optimal setting of initial weights. The current best practice seems to be Xavier initialization (see here for an informal explanation. Use Google or look at the original paper if you want the details)

Choice of loss function

In the example above I just chose to minimise squared loss but this is not the only choice. I don't want to go into further discussion of other options here, but I should be clear that squared loss is not the only (or even the best) choice.

Regularisation, Test/Training set separation, yadda yadda

In the interests of keeping the example simple, I left out a lot of stuff. It doesnt really make sense, for example, just to train on a training data set. In reality we want to use a test set to tell us when we are overfitting our data, and maybe even a validation set to see how well we will do on 'unseen' data. I've deliberately left out all this stuff to keep things simple in the example. Similarly, it sometimes makes sense to include some regularisation in the loss function -- maybe penalising larger weights. This helps prevent overfitting.

Adding noise to input

Another idea to improve the auto-encoder and prevent it overfitting is to add noise to the input, but try to get the network to learn/predict the input without the noise. The justification here is that we want the network to learn the 'essential' features of the input, and by adding noise to the input (but not the target/output) we are encouraging it to do this, rather than just memorize the input.

Getting started with deep learning in python...installation of keras with tensorflow backend

I've put these instructions down for installing keras to work with tensorflow in case they are helpful for others.

I'm going to assume that you are already on top of how to install python packages, or else have the excellent Anaconda (from continuum analytics) installed, which makes it easy via it's conda package manager (You can get it at https://www.continuum.io/).

If you have Anaconda installed, you should use its built-in package manager to install additional packages needed by keras:

conda install pyyaml
conda install h5py

You probably also want to make sure you are up to date with scipy, numpy, six:

conda upgrade six
conda upgrade scipy
conda upgrade numpy

Now, keras needs a back-end of either Theano or Tensorflow. I started using Tensorflow (for no particularly good reason, that's just what I picked). I just followed the Tensorflow setup instructions for Mac, which essentially consisted of typing the following at the command prompt:

pip install --upgrade https://storage.googleapis.com/tensorflow/mac/tensorflow-0.6.0-py2-none-any.whl

Now, finally, you want to install keras, which is the python library that sits atop Tensorflow/Theano and makes it easier to specify and train a neural network. You can clone the git repository which lives here. I just downloaded the zip file (from the same location), and unzipped it.

The keras website says to cd to the keras directory and run sudo python setup.py install

However, this did not work for me, because keras attempts to install Theano, and I had problems with that part of the install. Theano is listed as a dependency I suspect for historical reasons, even if you intend to use keras with a tensorflow backend. So, I just edited setup.py and removed theano as a dependency. By the time you read this, this issue will probably be fixed, and you wont need to edit setup.py.

Lastly, before you run keras, you need to create a .keras directory in your home folder, and in that directory, create a file called keras.json, and put the following in it: {"epsilon": 1e-07, "floatx": "float32", "backend": "tensorflow"}

Now, you should be able to run keras. At least this seems to have worked for me.

Public accountability gone in motorway funding

On the weekend the state government announced the "NorthConnex": another tolled motorway for Sydney. This one will cost $3 billion. That's small relative to the $11 billion WestConnex, but hardly small change. Some people think motorways are a bad idea because they induce more traffic, others think it is folly to invest billions of dollars in roadways when car travel is plateauing or falling in developed countries. Personally, I'm not opposed to motorways provided that:

1. They have a benefit/cost ratio significantly above one
2. The assumptions underlying the cost/benefit study are reasonable (there are many ways to game these things)

And this brings me to my main issue with these motorway proposals: there is no publicly available cost/benefit study for either of them. In a sane world, before we agreed to investing billions of dollars in infrastructure, the case for doing so would be publicly laid out. But it isn't, and the politicians are doing society two great disservices by their failure to do this: firstly, they may be wasting our money on ill-conceived projects; and secondly, they are entrenching a culture of opaqueness in the handling of public funds. This second point is especially important in NSW given our terrible record on corruption: if politicians can spend billions of dollars without the accountability of a public cost/benefit study, there will always be the potential for politicians to approve dud projects for a kickback, or a position on the board of a road construction or motorway operating company on retirement. Just so I am clear, I am not claiming that the WestConnex and NorthConnex motorways are going ahead due to corruption, I am simply pointing out that there are no publicly available cost/benefit studies, and we are being asked to trust the politicians and motorway companies without proof. This is no way to run major infrastructure investment.

My second peeve about these projects is that both of them are relying, to a significant degree, on tolls to fund them. This is a bad idea because funding a motorway this way lowers the benefit/cost ratio. This is easy to show mathematically, but is intuitively explained as follows:

• The motorway operator has an incentive to maximize cost-recovery/revenue, and the toll that does this is not the same as the toll that maximizes public benefit.
• Tolling individual links in a road network causes perverse outcomes: witness the Cross-City tunnel, which is tolled, and under-utilised, while the equivalent untolled surface route is still congested.

Marginal cost of congestion in Sydney is around $0.7/km

Post Synopsis: Every kilometre you drive in Sydney during peak hour, you are costing other motorists $0.66.... so if you drive 20km to work, other people are paying for your trip to the tune of about $13. If that sounds OK to you, remember that when you are stuck in traffic, you are paying (in time, extra fuel, etc) for someone else's decision to drive during peak hour. The best way out of this mess is congestion charging.

Congestion sucks

We all know that congestion is a drag in Sydney, and in most other large developed and developing cities (for fun, check out this congestion in Sao Paulo). A couple of studies have had a stab at estimating the total cost of congestion in Sydney. They found the following:

• A CIE report estimated congestion costs of ~ $0.29 per vehicle kilometre in 2005. These are total costs (not avoidable/dead-weight)
• BITRE estimate congestion costs of ~ $0.10 per vehicle kilometre in a 2007 working paper (a really good paper, btw). These are avoidable/dead-weight costs[*] -- total costs are around double that, so not so far from the CIE study.

I've reported these results (and others, below) in costs per vehicle-km, because thats the most common tangible measure, but if you want the total annual cost, you can multiply the above by the number of vehicle-km in Sydney (42.4 billion in 2005, say), and get total annual costs of $8-12 billion/year (depending on whether you use BITRE or CIE total costs).

Why average congestion costs are not so useful

Having just told you about a couple of prior studies that give average congestion costs, I have to break to it you that average costs aren't really that useful. Why? Two reasons:

• Firstly, there is that old problem with averages: your average temperature can be fine if you have your head in the oven and your feet in the freezer....  If you are driving during a relatively quiet period, when congestion is low or non-existent, then you aren't imposing costs on others. If, on the other hand, you are driving during peak hour, then you are going to be imposing very high costs. 
• Secondly, congestion is 'caused' by those last few thousand cars that are on the network. People refer to this effect as the 'school holiday effect', because you only need to remove a small amount of traffic for congestion to improve. This is telling us that the marginal cost of each additional vehicle is much higher than the average cost across all vehicles.

So what really matters, from a policy point of view, is the marginal cost of congestion. That is, what is the cost imposed by each additional driver? Since congestion costs are high during the peak, in this blog post I'm going to work out a rough estimate of the marginal congestion cost of each additional vehicle during peak hour in Sydney. 

So what are the marginal costs?

We are interested in the costs imposed on other motorists by each additional driver at peak time. We can work this out, to a reasonable approximation, as follows:

The average traffic speed during the inter-peak period is 42.9km/h [**]. The average traffic speed during the peak period is 30.0 km/h [***], so this means that traffic is 43% faster in the inter-peak period than during the (congested) peak period.

What causes this? Extra traffic on the roads, of course. Looking at the Figure below, we can see that during peak times, the volume of traffic on Sydney roads is ~ 270,000 vehicles during peak times, but only ~160,000 during inter-peak.

Source: BTS HTS 2011/12 summary report.

So, an extra 110,000 vehicles on the road network increases trip times by 43%. On a per-vehicle basis, this means each additional vehicle slows down all other traffic by 0.000325%. So if the vehicle is on the road for a minute, it slows down all other traffic by 0.000325% for that minute. Since there are 270,000 vehicles on the road at peak time, this means it costs those other motorists 0.88 minutes of delay. Valuing that at $14/hour (a pretty standard choice for peak travel), that's $0.205 of delay costs per minute of peak travel. If we add in other costs (extra running costs, pollution, etc), as in the BITRE report, this comes up to $0.33/minute. Average peak travel speed is 30 km/h, so converting this to a per-km costs, we get $0.66/km.

Note that this figure (of $0.66/km) is an average across Sydney for the peak period, and would be much higher in specific (congested) parts of the road network (and lower elsewhere). I haven't done the numbers, but my guess is that it would be a few times higher on congested links, which implies the marginal costs on congested roads might be in the vicinity of $2/km in the peak period.

[*] In the economics of congestion, there is an 'optimal' level of congestion: each additional traveller derives a benefit from making a trip, but there is a cost to other motorists due to congestion, and the optimal level of traffic is when the benefit to the additional traveller is equal to the costs imposed on existing travellers. Avoidable/dead-weight-loss congestion costs are the net costs taking into account the benefits of additional travel as well as the costs. Total costs are usually just the costs of delay relative to free-flow speeds.

[**] I've calculated average speed by using inter-peak travel times and distances from the NSW BTS's Strategic Travel Model. I've calculated the average speed across all trips in the ABS's journey-to-work data-set, using these travel times/distances.

[***] calculated as in [**], but I use am-peak travel times/distances.

Computing polygon centroids for non-contiguous polygons in R

I don't like doing spatial analysis in a graphical GIS program like Arc or QGis. I find them fine for looking at data, but not for analysis. For that, I'd rather use a scripting language, and I am NOT going to use Arc's python scripting support (my experience with arcgisscripting is not good).

So, I'm quite excited about the growing support for spatial data in Python and R (Yes, I probably need to get out more). I would love to get to the point where I could do all my analysis in Python or R and dump out the results for viewing in a GIS program.

So, in order to learn about the new R packages for spatial data, I set myself the task of calculating the centroids of non-contiguous polygons. So, more concretely, I mean:

• I have a layer of polygons, and each polygon has a name/identifier (a postcode, zone number, etc)
• I would like to calculate the centroid of each name/identifier. Note that there may be more than one polygon with that identifier, and in that case I want the centroid of all the polygons with that name/identifier. (Or, if you prefer, the centre of mass of all the polygons with the same identifier).

Note that the 'centroid' I am calculating here for each name/identifier may not be inside a polygon with that name/identifier. 

And here is the R code to do it (NB: In this data-set, the 'id' field for each polygon is 'tz06'):

#import everything you need in 1 package

library(GISTools)

#read in your shapefile

p = readShapePoly("GMA_polygons_wgs84")

#calculate the centroids of each polygon

trueCentroids = gCentroid(p, byid=TRUE)

#I like to whack the centroids in as extra columns

p$lat = coordinates(trueCentroids)[,2]

p$lon = coordinates(trueCentroids)[,1]

#The 'average' lat and lon is what I want. It is the centre 

#of mass of all polygons with that id. 

#Of course, if there is only a single polygon with that id,

#the average lat/lon is just the regular centroid, so I 

#initialize my centroid to that first

p$avglat = p$lat

p$avglon = p$lon

#Now go through and calculate the 'average' centroid

#for all the tricky areas that are made up of more than

#1 polygon

for(i in 1:length(p$tz06)) 

{

    #find all the polygons with the same id

    matches = p$tz06 == p$tz06[i]

        

    #if there is more than one such polygon

    if(sum(matches) > 1)

    {

        #we weight all of the centroids by area,

        #because we want the centre of mass

        areas = matches*p$Area_sqm

        weights = areas / sum(areas)

        p$avglat[i] = sum(p$lat*weights)

        p$avglon[i] = sum(p$lon*weights)

    }

}

I'm sure this is not the nicest way of doing this, so if you know a cleaner way (hell, there may even be a built-in function for this already!), please let me know. Also, I'm hopeless with regards to functional programming (my brain thinks procedurally), so I've dont the above with a loop, but if you can point out to me a nice functional way of achieving the above with map/filter/fold/etc, I'd be keen to hear about that too.

The disaster that is Google Drive

My wife is a smart lady. Not particularly interested in technology, but no luddite either. She recently signed up with Google Drive so that pictures and videos of our kids could be backed up online.

Big mistake.

Files mysteriously go missing from our home computer's LOCAL folders (i.e. the folder of images and videos she wants to sync). In fact many of the local folders were empty! Some have strange duplicates (for example, a "Pictures" and a "Pictures(1)" folder). This was quite a shock when she first discovered it -- "Google Drive is deleting the precious files I want to backup!!!". Looking into it, we discovered that other people have had the same problem.

We were able to find some of our 'missing' files in the "Trash" on our local PC. This caused us even more alarm, because we still lost a lot of really precious stuff from our local computer. So what this meant for us is that while we were able to go to Trash and recover some of the files, we lost many of them. My wife also assured me she wasn't doing anything strange like moving stuff around and deleting stuff while trying to back it up.

This seems really bizarre behaviour for a cloud storage solution. It should never delete stuff from your local computer for you. Or perhaps there are a few cases where it should (if you have manually deleted those items on another computer that syncs with your drive account), but if it does, it should be really careful about it. My wife didn't do anything to suggest to Google Drive that it should delete those items.


So, in panic we shut off Google Drive on our home computer and went to the web interface to see what was there. Alarmingly, all the folders we looked at were empty! Real panic stations now... but looking at our usage, we found that we were using ~20GB of Drive storage. Digging a little further, I discovered that I could click on the (hard to find) "All Items" link on the web interface and we could see everything that was uploaded (actually I think this means that it was moved to the online Google Drive Trash, but online Trash fortunately has to be emptied manually). I manually downloaded these (the interface to do this is terrible, as an added bonus, and you can only download in 2GB chunks). End result is that I think I've recovered our stuff. But the first thing I will do when I've confirmed that I've recovered everything is delete our Google Drive account and sign up with Box or DropBox or someone (anyone!) else.


Note that I am a big user of Google products, and I'm generally pretty disposed to liking Google: my work uses gmail for corporate email, I use it for my personal email, share videos via youtube, and so on. Heck this blog is run by Google. But Drive has been a massive headache for us.

Thinking of using Google Drive? Just dont do it.

How to make high-speed rail stack-up

On the face of it, the latest of a long line of studies into fast rail appears to show, yet again, that it is going to be hard to justify the huge up-front capital costs. An estimated cost of $80 billion (give or take $20-odd billion, but who cares about loose change?), which is not going to be recouped in any significant way from fares. So lets not bother, right?

But hang on, aren't we missing something here? Well, only the main cost-recovery avenue -- value (re)capture. The idea is hardly new, but for some reason I can't fathom it is rarely discussed in detail in these fast rail feasibility studies. For those unfamiliar with the idea, it's quite simple.  The government is going to shell out huge buckets of money building a high-speed rail line. Anyone owning property near the rail stations will immediately see the value of that property increase. Since these property owners haven't done anything themselves, and are just benefiting from the government's spending, it's fair for the government to recoup the increase. At this point people often start quibbling about difficulties in estimating the benefits, or object to trying to claw back money from old ladies who just happen to own a house near one of the stations. None of that is relevant except in political terms, and for the moment lets stay away from the politics: we can come back to that later.

So how important is value-capture, relative to fare-box revenue? Very. Lets do our own back-of-the envelope calculations and see. First up, lets simplify things just to make the case for value recapture clear. Let's just take the Sydney-Newcastle-Canberra section of the track. The best-guess cost for this in the most recent study is $29.3 billion; let's just call it $30 billion.

Total cost: $30 billion

A further simplification will help bring out my point, and allow us to avoid complicated calculations that just obscure my main point. Let's say that instead of linking Sydney, Canberra and Newcastle, we just build south-west from Sydney out to empty sheep paddocks near Golbourne, and North to some easily developable land near the Central Coast or Newcastle. So, we have two lines, with 4 stops along each line (as in the most recent study), and a train commute of 40 minutes or under for each station.

Now, suppose we develop a small 2km radius 'town' around each station. Since we've built out to sheep paddocks this should be no problem. A 2km radius town with 8 lots per hectare will give us 10000 lots in each town. How much would each lot be worth? With a sub-40minute commute to Sydney, I'd say at least three or four hundred thousand. But lets say after servicing and other costs, the net increase is $200k for each lot. So, 8 stations, 10000 lots per station, $200k per lot, gives us $16 billion dollars:

Land-capture capital recovery: conservatively, $16 billion

So with some pretty conservative assumptions we have recouped half the capital cost. Getting a little more creative should allow even greater cost recovery. Developing a 3km^2 town at each station, or developing at higher densities, are easy adjustments that would increase revenue and get you to near full cost-recovery without any fare-box revenue.

So what about fare-box revenue? In our 8 towns, with 10000 dwellings, lets say there are 10000 workers, 30% of which commute using the rail line, with the other 70% commuting by car or some other mode. So, 3000 workers in each town commuting each way 200 days a year gives us 9.6 million trips. Call it 10 million. Add another 10 million non-work trips. Say each trip costs $10. So, that's $200 million a year, or a 0.6% p.a. gross return on $30 billion invested, before we take into account operating costs. Fiddle with the above numbers any way you like, it wont change the fact that fare-box revenue is almost small enough to forget about altogether, and certainly much less important than value-capture as a funding mechanism.

Net fare-box revenue: approximately nil.

Bottom-line: if you want to fund high-speed rail, heck, in fact rail of any kind, you need to put value recapture front and centre of your approach to cost-recovery. Why isn't it already? Because property owners would prefer to pocket the benefits from the rail-line themselves, and many will fight tooth-and-nail to do so. Building the line to some empty sheep paddock mostly avoids this, because there are only a handful of large land-owners likely to be involved, and that's why I took that approach in the above example, but in reality there would be a lot of opposition to such a move ("It's crazy to build a rail line out to nowhere! We need to strengthen the centres we already have!"). And  if you do build to existing centres, it is politically difficult to claw back the increased land values (not that any Australian government have really tried very hard....). In fact, building the line out to empty paddocks is probably not the best plan, economically, anyway, provided that you have a mechanism to recoup the increased value from existing land-holders.

And the reason land-value capture works so well (if you can do it, politically), is because our dysfunctional planning systems have made supply of new dwellings, at low cost, so difficult. So there is a lot of suppressed demand for accessible locations. A rail line is, essentially, a way of 'manufacturing' new accessible locations. Releasing land on the urban fringe is not the same thing -- a point missed by many who claim that the fix to our housing supply problems is just to release more land. This would help a bit, but the real demand is not for land, but for dwellings with good accessibility (to jobs, and other things), and accessibility can only really be supplied by transport infrastructure.

As a final aside, some economists have pointed out that because the benefits of rail systems get capitalized into a relatively small number of land-holders, this makes rail distributively unfair, as costs are shared and benefits and concentrated. The benefits of a motorway, on the other-hand, are much more widely spread. This is a valid point, and can't be addressed without re-thinking value capture. Just to be clear: I'm not arguing for motorways here, just making it clear that making rail work requires us to have a discussion about the fairness of recouping private gains from public infrastructure investment.

Connecting to mysql from python the odbc way from mac

This post just describes how to set up python to talk to a mysql database on a mac. Written for reference, and for co-workers who will shortly need to do something like this.

I wanted to be able, in a python script, to talk to a remote MySQL database from my mac. As I went about this I found the documentation a little scattered -- I couldnt find a single page giving the instructions, and so I've just written up the steps here.

I assume that you have python installed already.

Step 1: Install pyodbc

First step is to install pyodbc if you havent got it installed already. I downloaded pyodbc from http://code.google.com/p/pyodbc/downloads/list (I downloaded the source distribution because there were no binaries for mac. The source distribution is the one that ends in .zip)

After unzipping this, going into the directory that is created by the unzip and typing python setup.py install at the command prompt should install everything without you needing to do anything else.

Step 2: Check that you have an odbc manager installed

pyodbc wont do all the work for you -- it relies on an odbc manager to ferry things back and forth to databases. On a Mac you should have an odbc manager installed already by default. My machine has one called iodbc: if I type iodbctest at the command prompt* I get the following:

iODBC Demonstration program
This program shows an interactive SQL processor
Driver Manager: 03.52.0607.1008

Enter ODBC connect string (? shows list):

If I then type '?' as suggested, I get an empty list... which explains that there are no drivers installed. What is a driver? Its a database-specific connector that lets the odbc manager talk to a particular database (MySQL, Postgresql, SQLServer, etc). So the last link required in the communication chain is a mysql driver.

* Apparently you can get to iodbc through Applications->Utilities as well if you want to look for it in the Mac GUI.

Step 3: Install your odbc drivers

OK, so we now need to install specific odbc drivers for each database we want to connect to. In this case, I'm just going to install MySQL. I grab the odbc driver from http://dev.mysql.com/downloads/connector/odbc (picking the Mac .dmg file just to make the installation easy). A few clicks on the .dmg file and the driver is installed.

I can check that it is installed by running iodbctest again from the command prompt. This time when I type '?' to see what drivers are installed I get:

Enter ODBC connect string (? shows list): ?

DSN | Driver
------------------------------------------------------------------------------
myodbc | MySQL ODBC 5.1 Driver

OK, so that completes the pipeline. Now in python I can use pyodbc, which talks to iodbc (or whatever ODBC manager you have installed), which uses the mysql odbc driver to talk to a MySQL database.

So lets do that and check that it works

Step 4: connect to your database from python

OK, so now we should be able to connect to your database in python.

But one last thing you need to do is to work out what connection string you need to connect to your database. A helpful resource with example connection strings can be found at http://www.connectionstrings.com/. Here is what I did in python to connect to my database:

>>> import pyodbc
>>> connstr="Driver={MySQL ODBC 5.1 Driver};Server=your.server.;Port=3306;Database=databasename;User=username; Password=nottelling;Option=3;"
>>> conn=pyodbc.connect(connstr)
>>>

and from here you can just play with the conn object as you like in python! So, for example, to look at a table:

>>> res = conn.execute("select * from sometable")
>>>for r in res:
>>>    print r.next()

Rental yields will increase, but not the way they'd like....

In response a growing chorus of analysts pointing out that Australian property seems to have both limited prospects for capital growth and low rental yields, many property market boosters respond that yields will increase (even the boosters have given up spruiking capital growth as a reason to buy property -- that boat has sailed).

Now, for once I am in agreement with the property markets boosters -- yields must increase. But there are two ways for yields on residential property to increase. Increasing rents is one (and clearly the one the boosters have in mind). But falling prices will do the job just nicely also. Just saying.

newspapers snoozpapers

On a recent flight back from Perth, I found myself growing increasingly exasperated at the poor quality of the journalism in the two newspapers available -- the Fin Review and the Australian.

For some time I've been hoping, against all evidence to the contrary, that the Fin would not follow the Sydney Morning Herald's slide into mindless he-said she-said reporting. But reading the Fin during the last week has not been encouraging. Why read the Fin if not for some actual analysis? But what we get now on almost every issue are the verbatim opinions of pundits -- financial, political, social. Knowing what Bill Gross' views are on inflation and government indebtedness, while interesting, is complementary to, not a replacement for, detailed reporting.

The Saturday article on 'bond vigilantes' was especially unsatisfying -- with extensive quotes from 'senior bond traders' together with repetition of all the tired old news-bites about looming US inflation and sovereign debt concerns. Did the article contain any references to actual evidence or analysis enabling readers to evaluate such claims? US core inflation figures? Long term bond yields consistent with either rampant inflation or sovereign default? No, just appeals to financial market authority-figures.

This low-quality style of reporting has been prominent in real-estate journalism for a long time, where puff pieces quoting real estate agent and industry mouthpieces have forced those interested in any thoughtful analysis to other sources. (Leith van Onselin, for example, has posted manytimes on the issue -- such as here, here, and here).

As for the Australian, the problems there are quite different. Not being a regular reader, I was unprepared for the breadth and depth of the editorial lunacy, from Greg Sheridan's column railing against the economy wrecking carbon tax to Chris Merritt's nanny-state rant. In political and news sections, Peter van Onselin's column offered the only relief (Peter is Leith's brother -- there must be something in the genes). In Business, Glenda Korporaal's piece was also thoughtful and balanced enough to be out-of-place.

Not content with its plentiful home-grown lunacy, the Australian even found space to import some, with William McGurn's WSJ piece praising Paul Ryan's fairytale deficit reduction plan -- entirely neglecting to mention that the plan relies entirely on heroic assumptions (2.8% unemployment in 2021, anyone?).

I'm too young to become a grumpy old man about these things, but they're forcing my hand dammit!

Just a reminder to myself on converting between coordinate systems in ArcGIS.

I often want to convert data from one coordinate system into another, and I frequently forget how to do this. This is just a post to myself so that I remember how to do this..

Below is just a cached copy of this page, replicated here in case the original page goes away.....

---

Converting your vector files into different Coordinate Systems Using ArcGIS 8.2

EXAMPLE: converting NAD83 MTM ZONE8 to WGS84 UTM ZONE18

1. open ArcToolBox
2. choose Data Management Tools → Projections → Project Wizard (shapefiles, geodatabase)

3. select the shapefile to be converted, then its original coordinate system will be list in the dialog. Click Next.

4. define output file name and stored location; click NEXT.

5. giving desired coordinate system; click Select Coordinates System

6. click Select to define a new coordinate system

7. choosing Geographic (lat/long) or Projected (x/y) coordinate system.

For WGS84 ZONE18N, choose Projected coordinate system → UTM → WGS84 → WGS84 UTM ZONE18N. Click NEXT.

8. Detail information of WGS84 UTM ZONE 18N will be list in the dialog. Click NEXT

9. A Geographic transformation (optional) will be required since (in this case) NAD27 and WGS84 are different datum. A list of option will be offered in the function; choose the newest transformation (or fit the case). Finally, your output file should be projected in the new coordinate system.

Schumacher's foresight and the folly of complex forecasting models

I've been reading Schumacher's "Small is Beautiful" recently -- quite an old book now but still worth a read. He writes very clearly, and the subject matter of the book is mostly still relevant, and the key ideas still thought-provoking. He shows uncommon foresight on a few topics. For example, he talks about the near certainly of a future where oil is scarce & expensive, and he was doing so before even the 1970's oil shocks. He also talks about the folly of eroding natural capital long before it was common to do this. But perhaps most interestingly, he identifies the process where a void in values (created by the decline of religion) was being filled by economics, and despite being an economist, he thought this dangerous. He thought it far safer to rely on traditional cultural & religious value systems than to subject everything to economic cost benefit analysis. I think this last point is especially interesting because, while most economists will probably baulk at it, recent events lend weight to Schumacher's view. For example, look at all the economic debate and analysis that has gone into climate change, and look where it has got us. Wouldn't we have done better developing a value system where caring for our planet was more than just an economic consideration, but something sacred?

Schumacher wasnt a hippy or anything either, nor was he a fuzzy-minded economist who hated the reduction of everything to numbers because he himself couldnt do any maths. In fact, he was a competent mathematician. This makes some of what he has to say about model-building in the later chapters of his book particularly interesting. Being a bit of a data-cruncher and model builder myself, I liked the quotes in particular, which are right on the money, and which I think should be compulsory reading for any model builder:

"Crude methods of forecasting ... are not likely to lead into the errors of spurious verisimilitude and spurious detailing -- the two greatest vices of the statistician. Once you have a formula and an electronic computer, there is an awful temptation to squeeze the lemon until it is dry and present a picture of the future which through its very precision and verisimilitude carries conviction. Yet a man who uses an imaginary map, thinking it a true one, is likely to be worse off than someone with no map at all; for he will fail to inquire wherever he can, to observe every detail on his way, and to search continuously with all his senses and all his intelligence for indications of where he should go.

The person who makes the forecasts may still have a precise appreciation of the assumptions on which they are based. But the person who uses the forecasts may have no idea at all that the whole edifice, as is often the case, stands and falls with one single, unverifiable assumption. He is impressed by the thoroughness of the job done, by the fact that everything seems to 'add up;, and so forth. If the forecasts were presented quite artlessly, as it were, on the back of an envelope, he would have a much better chance of appreciating their tenuous character and the fact that, forecasts or no forecasts, someone has to take an entrepreneurial decision about the unknown future."

I wish I could say I had never fallen for the temptation to go for a complex model with 'rich' outputs, when it would have been essentially as accurate, and less misleading, to present a few back of the envelope projections.....