Friday, February 14, 2014

capturing stdout in a python child process

Lately at work I've been putting together a web management console that handles execution of some of our long running Gurobi models. This makes a lot of sense to do since the models run on Gurobi Cloud for hours or even days before achieving their desired optimality gaps. Gurobi Cloud requires a persistent connection, so we use a server to handle that instead of solving models from our desktop machines. In a sense, this is a bit of a throwback to the time sharing systems from before the PC era.

As a result, I've had to update my dormant (and, to be honest, never particularly acute) administration and systems programming skills. There are always some interesting gotchas that come up in this sort of work. This time I ran into the following conundrum: I have a cron job that is responsible for constructing model instances in a child process, handing them off to Gurobi Cloud, and saving their results. Its output is redirected to a file so I can tail it in real time and store it for later. For the web console, I also want to store the output to a database so I can easily display it. And I'd like to do this in as generic a manner as possible, so I can add other types of models (custom heuristics, for instance) and also capture anything they write to standard output and standard error.

There are a number of threads available on capturing standard output. Some of the advice I found recommended using os.dup2 to clone the stdout file descriptor into another file. The problem with this approach is it requires having an open file, which is inconvenient, and it ends up closing the existing stdout. So the output does get redirected, but then it disappears from my log. Also, I could only seem to capture it once. I wanted to do this repeatedly.

Other advice got me closer to the mark: Replace sys.stdout with an instance of StringIO, then just read its contents later. This allowed me to do redirections in serial, but once again meant that standard output was essentially shut off.

The final solution was to subclass StringIO and override its write method. Every time write is called, it writes to both the in-memory string and a canonical reference to stdout. Example code is given below that captures stdout and stderr from a child process while still letting their messages be visible.

from StringIO import StringIO
from multiprocessing import Process, Queue
import sys

# Proxy classes that both write to and capture strings
# sent to stdout or stderr, respectively.
class TeeOut(StringIO):
    def write(self, s):
        StringIO.write(self, s)
        sys.__stdout__.write(s)

class TeeErr(StringIO):
    def write(self, s):
        StringIO.write(self, s)
        sys.__stderr__.write(s)

class Parent(object):
    def run(self):
        # Start child process.
        queue = Queue()
        child = Process(target=self.spawn, args=(queue,))
        child.start()
        child.join()

        # Do someting with out and err...
        out = queue.get()
        err = queue.get()

    def spawn(self, queue):
        # Save everything that would otherwise go to stdout.
        out = TeeOut()
        sys.stdout = out

        err = TeeErr()
        sys.stderr = err

        try:
            # Do something...
            print 'out 1'
            print 'out 2'
            sys.stderr.write('error 1\n')
            print 'out 3'
            sys.stderr.write('error 2\n')

        finally:
            # Restore stdout, stderr and send their contents.
            sys.stdout = sys.__stdout__
            sys.stderr = sys.__stderr__
            queue.put(out.getvalue())
            queue.put(err.getvalue())

if __name__ == '__main__':
    Parent().run()

A few minor points:

  • StringIO can't be subclassed like object-based types, thus lines 9 and 14 can't use super.
  • The sys module keeps two references to stdout: sys.stdout and sys.__stdout__. We use the latter to recover the former after processing.
  • The contents of standard out and standard error are transmitted from the child process to its parent at the end of processing. A smarter implementation might transfer chunks of them periodically, but I'm not expecting them to be large.
  • The transfer happens in a finally block. Thus any run time exception in the child will not interfere with logging.

Monday, February 3, 2014

chebyshev centers of polygons with gurobipy

A common problem in handling geometric data is determining the center of a given polygon. This is not quite so easy as it sounds as there is not a single definition of center that makes sense in all cases. For instance, sometimes computing the center of a polygon's bounding box may be sufficient. In some instances this may give a point on an edge (consider a right triangle). If the given polygon is non-convex, that point may not even be inside or on its boundary.

This post looks at computing Chebyshev centers for arbitrary convex polygons. We employ essentially the same model as in Boyd & Vandenberghe's Convex Optimization text, but using Gurobi instead of CVXOPT.

Consider a polygon defined by the intersection of a finite number of half-spaces, $Au \le b$. We assume we are given the set of vertices, $V$, in clockwise order around the polygon. $E$ is the set of edges connecting these vertices. Each edge in $E$ defines a boundary of the half-space $a_i'u \le b_i$.

(x, y)
$V = \{(1,1), (2,5), (5,4), (6,2), (4,1)\}$
$E = \{((1,1),(2,5)), ((2,5),(5,4)), ((5,4),(6,2)), ((6,2),(4,1)), ((4,1),(1,1))\}$

The Chebyshev center of this polygon is the center point $(x, y)$ of the maximum radius inscribed circle. That is, if we can find the largest circle that will fit inside our polygon without going outside its boundary, its center is the point we are looking for. Our decision variables are the center $(x, y)$ and the maximum inscribed radius, $r$.

In order to do this, we consider the edges independently. The long line segment below shows an arbitrary edge, $a_i'u \le b_i$. The short line connected to it is orthogonal in the direction $a$. $(x, y)$ satisfies the inequality.

(x, y)

The shortest distance from $(x, y)$ will be in the direction of $a$. We'll call this distance $r$. If we were to move the edge so it had the same slope but went through $(x, y)$, its distance from $a_i'u = b_i$ would be $r||a_i||_2$. Thus we can add a constraint of the form $a_i'u + r||a_i||_2 \le b_i$ for each edge and maximize the value of $r$ as our objective function.

max $r$
subject to $(y_i-y_j)x + (x_j-x_i)y + r\sqrt{(x_j-x_i)^2 + (y_j-y_i)^2} \le (y_i-y_j)x_i + (x_j-x_i)y_i\, \forall\, ((x_i,y_i), (x_j,y_j)) \in E$

As this is linear, we can solve it using any LP solver. The following code does so with Gurobi.

#!/usr/bin/env python
from gurobipy import Model, GRB
from math import sqrt

vertices = [(1,1), (2,5), (5,4), (6,2), (4,1)]
edges = zip(vertices, vertices[1:] + [vertices[0]])

m = Model()
r = m.addVar()
x = m.addVar(lb=-GRB.INFINITY)
y = m.addVar(lb=-GRB.INFINITY)
m.update()

for (x1, y1), (x2, y2) in edges:
    dx = x2 - x1
    dy = y2 - y1
    m.addConstr((dx*y - dy*x) + (r * sqrt(dx**2 + dy**2)) <= dx*y1 - dy*x1)

m.setObjective(r, GRB.MAXIMIZE)
m.optimize()

print 'r = %.04f' % r.x
print '(x, y) = (%.04f, %.04f)' % (x.x, y.x)

The model output shows our center and its maximum inscribed radius.

r = 1.7466
(x, y) = (3.2370, 2.7466)
(x, y)

Question for the reader: in certain circumstances, such as rectangles, the Chebyshev center is ambiguous. How might one get around this ambiguity?

Monday, July 29, 2013

network splitting

Last week, Paul Rubin wrote an excellent post on Extracting a Connected Graph from an existing graph. Lately I've been performing related functions on data from OpenStreetMap, though without access to a solver. In my case I'm taking in arbitrary network data and splitting it into disconnected subnetworks. I thought it might be a good case study to show an algorithmic way doing this and some of the performance issues I ran into.

A small example can be seen below (click on the image to enlarge it). This shows a road network around the Las Vegas strip. There is one main (weakly) connected network in black. The roads highlighted in red are disconnected from the main network. We want code that will split these into connected subnetworks.

Say we have data that looks like the following. Instead of nodes, the numbers in quotes represent edges. Think of these as streets.

{
    "0": [1, 2, 3],
    "1": [9248, 9249, 9250],
    "2": [589, 9665, 9667],
    "3": [0, 5, 6],
    "4": [0, 5, 6],
    "5": [588],
    "6": [4, 8, 9],
    ...
}

Our basic strategy is the following:

  1. Start with every edge alone in its own subnetwork.
  2. For each connection, merge the networks of the source and destination edges.

#!/usr/bin/env python
import json
import sys
import time

class hset(set):
    '''A hashable set. Note that it only hashes by the pointer, and not by the elements.'''
    def __hash__(self):
        return hash(id(self))

    def __cmp__(self, other):
        return cmp(id(self), id(other))

if __name__ == '__main__':
    try:
        inputfile = sys.argv[1]
    except:
        print 'usage: %s network.json' % sys.argv[0]
        sys.exit()

    print time.asctime(), 'parsing json input'
    connections = json.load(open(inputfile))

    edge_to_net = {} # Edge ID -> set([edges that are in the same network])
    nets = set()     # Set of known networks

    print time.asctime(), 'detecting disconnected subgraphs'
    for i, (from_edge, to_set) in enumerate(connections.iteritems()):
        from_edge = int(from_edge)

        try:
            from_net = edge_to_net[from_edge]
        except KeyError:
            from_net = edge_to_net[from_edge] = hset([from_edge])
            nets.add(from_net)

        if not (i+1) % (25 * 1000):
            print time.asctime(), '%d edges processed / %d current subnets' % (i+1, len(nets))
        
        for to in to_set:
            try:
                to_net = edge_to_net[to]

                # If we get here, merge the to_net into the from_net.
                if to_net is not from_net:
                    to_net.update(from_net)
                    for e in from_net:
                        edge_to_net[e] = to_net
                    nets.remove(from_net)
                    from_net = to_net

            except KeyError:
                from_net.add(to)
                edge_to_net[to] = from_net

    print time.asctime(), len(nets), 'subnets found'

We run this against the network pictured above and it works reasonably quickly, finishing in about 7 seconds:

Mon Jul 29 12:22:38 2013 parsing json input
Mon Jul 29 12:22:38 2013 detecting disconnected subgraphs
Mon Jul 29 12:22:38 2013 25000 edges processed / 1970 current subnets
Mon Jul 29 12:22:44 2013 50000 edges processed / 124 current subnets
Mon Jul 29 12:22:45 2013 60 subnets found    

However, when run against a road network for an entire city, the process continues for several hours. What is the issue?

The inefficiency occurs from lines 46 to 50. In this we are frequently removing references to every element in a large set. Instead, it would be better to remove as few references as possible. Therefore, instead of merging from_net into to_net, we will determine which network is the smaller of the two and marge that one into the larger one. Note that this does not necessarily change the worst case time complexity of the algorithm, but it should make the code fast enough to be useful. The new version appears below.

#!/usr/bin/env python
import json
import sys
import time

class hset(set):
    '''A hashable set. Note that it only hashes by the pointer, and not by the elements.'''
    def __hash__(self):
        return hash(id(self))

    def __cmp__(self, other):
        return cmp(id(self), id(other))

if __name__ == '__main__':
    try:
        inputfile = sys.argv[1]
    except:
        print 'usage: %s network.json' % sys.argv[0]
        sys.exit()

    print time.asctime(), 'parsing json input'
    connections = json.load(open(inputfile))

    edge_to_net = {} # Edge ID -> set([edges that are in the same network])
    nets = set()     # Set of known networks

    print time.asctime(), 'detecting disconnected subgraphs'
    for i, (from_edge, to_set) in enumerate(connections.iteritems()):
        from_edge = int(from_edge)

        try:
            from_net = edge_to_net[from_edge]
        except KeyError:
            from_net = edge_to_net[from_edge] = hset([from_edge])
            nets.add(from_net)

        if not (i+1) % (25 * 1000):
            print time.asctime(), '%d edges processed / %d current subnets' % (i+1, len(nets))
        
        for to in to_set:
            try:
                to_net = edge_to_net[to]

                # If we get here, merge the to_net into the from_net.
                if to_net is not from_net:
                    # Update references to and remove the smaller set for speed.
                    if len(to_net) < len(from_net):
                        smaller, larger = to_net, from_net
                    else:
                        smaller, larger = from_net, to_net

                    larger.update(smaller)
                    for e in smaller:
                        edge_to_net[e] = larger
                    nets.remove(smaller)
                    edge_to_net[to] = larger
                    from_net = larger

            except KeyError:
                from_net.add(to)
                edge_to_net[to] = from_net

    print time.asctime(), len(nets), 'subnets found'

Indeed, this is significantly faster. And on very large networks it runs in minutes instead of hours or days. On the small test case used for this post, it runs in under a second. While this could probably be done faster, that's actually good enough for right now.

Mon Jul 29 12:39:55 2013 parsing json input
Mon Jul 29 12:39:55 2013 detecting disconnected subgraphs
Mon Jul 29 12:39:55 2013 25000 edges processed / 1970 current subnets
Mon Jul 29 12:39:55 2013 50000 edges processed / 124 current subnets
Mon Jul 29 12:39:55 2013 60 subnets found

Wednesday, July 17, 2013

stack in go

I've been using an increasing amount of Go lately. One of my favorite features is its ability to add methods to any existing type. This is possible for primitive types by simply providing an alias for them, similar to using typedef in C.

Recently, I found I needed a stack container. Nothing complicated, merely an implementation of push, pop, and empty methods. While this is easy enough to do using slices, wrapping the required code makes for a much nicer interface.

package main

import "fmt"

type Element float64
type Stack []Element

func (s *Stack) Push(e Element) {
    *s = append(*s, e)
}

func (s *Stack) Pop() Element {
    e := (*s)[len(*s)-1]
    *s = (*s)[:len(*s)-1]
    return e
}

func (s *Stack) Empty() bool {
    return len(*s) < 1
}

func main() {
    s := Stack([]Element{})

    for i := 0; i < 10; i++ {
        s.Push(Element(i))
    }

    for !s.Empty() {
        fmt.Println("pop:", s.Pop())
    }
}

Since there aren't real classes in Go, or inheritance, this doesn't do anything scary like monkey patching. It's just a bit of syntactic sugar to make our lives better. For instance, the Push method could just as easily look like this to the compiler:

func Push(s *Stack, e Element) {
    *s = append(*s, e)
}

Monday, October 8, 2012

web slides with python & math using deck.js, syntaxhighlighter, and mathjax

In putting together some slides for class that include both Python code and math models, I realized that it was high time to update the technology stack I'm using. If anyone else finds this useful, please feel free to copy my starting template and use or improve upon it.

I'm hosting my slides directly out of Dropbox since I can do so without uploading them to a web server. It also means I can edit and view them locally without have to go through any sort of manual syncing process. This is possibly the simplest way in existence to share a page on the internet.

Steps:

  1. Create the directories /public/, /public/www/, /public/www/slides/, and /public/www/slides/boilerplate/ in your Dropbox folder.
  2. Download the latest version of deck.js and decompress it into /public/www/deck.js/.
  3. Download SyntaxHighlighter and decompress it into /public/www/syntaxhighlighter/.
  4. Download MathJax and decompress it into /public/www/mathjax/.
  5. Save this template into /public/www/slides/boilerplate/boilerplate.html.

Now to create a new slide deck, copy the boilerplate directory to something new in /public/www/slides/. You should be able to begin with something that looks like this. Search through its source for "slides start" to see how to create slides that show beautifully rendered Python and LaTeX.

Use the left and right arrow keys or the icons that appear when you hover over the left and right of the slides to navigate.

Monday, September 24, 2012

roll your own gurobi shell

Note: These instructions have been updated to use the setup.py from the Gurobi distribution.

I've recently started using Gurobi for a class. While I'm pleased with the level of Python support they provide, I'm not a fan of the Gurobi shell they ship it with.

  • The shell is Python, which is good. However, instead of bundling the Gurobi libraries as a Python module that I can install alongside my system's native Python, they ship their own Python binaries. If I want to use any third-party modules I have to make them available to Gurobi somehow.
  • They appear to have removed distutils from the Python they ship.
  • The Gurobi shell doesn't give me module docs, tab completion, or anything I'd get from either ipython or bpython. It's nothing more than the rudimentary cPython shell.

They do, however, provide a setup.py script, so we don't have to use gurobi.sh or the Python that's in the distribution. Here are simple instructions for using bpython as your shell. You should have Gurobi 5.0.1 installed and the GUROBI_HOME environment variable set. (I assume you are installing it into your home directory.)

$ mkdir -p $HOME/lib/python2.7/site-packages/
$ cd $GUROBI_HOME
$ PYTHONPATH=$HOME/lib/python2.7/site-packages/ python setup.py install --prefix=$HOME

Update your environment's PYTHONPATH and LD_LIBRARY_PATH:

$ echo "export PYTHONPATH=\$PYTHONPATH:\$HOME/lib/python2.7/site-packages/" >> $HOME/.bashrc
$ echo "export LD_LIBRARY_PATH=\$LD_LIBRARY_PATH:\$GUROBI_HOME/lib/" >> $HOME/.bashrc
$ source $HOME/.bashrc

Install bpython:

$ sudo apt-get install bpython

Now type bpython and you should have the lovely interface you see below, complete with tab completion and module docs. Isn't that better?

Saturday, September 22, 2012

lagrangian relaxation with gurobipy

Note: this is still a new technique to me, so please point out any errors and I will correct them in the post.

Correction: $b_{ij}$ is now $b_j$ in the second set of constraints.

We've been studying Lagrangian Relaxation (LR) in the combinatorial optimization course I'm taking this term, and I had some difficulty finding a simple example covering its application. In case anyone else finds it useful, I'm posting a Python version for solving the Generalized Assignment Problem (GAP). This won't discuss the theory of LR at all, just give example code using Gurobi 5.0.1.

The GAP as defined by Wolsey consists of a maximization problem subject to a set of set packing constraints followed by a set of knapsack constraints:

max $\sum_i \sum_j c_{ij} x_{ij}$
subject to $\sum_j x_{ij} \leq 1 \forall i$
$\sum_i a_{ij} x_{ij} \leq b_{ij} \forall j$
$x_{ij} \in \{0, 1\}$

An unadulterated version of this BIP using gurobipy might look like the following:

#!/usr/bin/env python

# This is the GAP per Wolsey, pg 182

from gurobipy import GRB, Model

model = Model('GAP per Wolsey')
model.modelSense = GRB.MAXIMIZE
model.setParam('OutputFlag', False) # turns off solver chatter

b = [15, 15, 15]
c = [
    [ 6, 10,  1],
    [12, 12,  5],
    [15,  4,  3],
    [10,  3,  9],
    [ 8,  9,  5]
]
a = [
    [ 5,  7,  2],
    [14,  8,  7],
    [10,  6, 12],
    [ 8,  4, 15],
    [ 6, 12,  5]
]

# x[i][j] = 1 if i is assigned to j
x = []
for i in range(len(c)):
    x_i = []
    for j in c[i]:
        x_i.append(model.addVar(vtype=GRB.BINARY))
    x.append(x_i)

# We have to update the model so it knows about new variables
model.update()

# sum j: x_ij <= 1 for all i
for x_i in x:
    model.addConstr(sum(x_i) <= 1)

# sum i: a_ij * x_ij <= b[j] for all j
for j in range(len(b)):
    model.addConstr(sum(a[i][j] * x[i][j] for i in range(len(x))) <= b[j])

# max sum i,j: c_ij * x_ij
model.setObjective(
    sum(
        sum(c_ij * x_ij for c_ij, x_ij in zip(c_i, x_i))
        for c_i, x_i in zip(c, x)
    )
)
model.optimize()

# Pull objective and variable values out of model
print 'objective =', model.objVal
print 'x = ['
for x_i in x:
    print '   ', [1 if x_ij.x >= 0.5 else 0 for x_ij in x_i]
print ']'

The solver quickly finds the following optimal solution of this toy problem:

objective = 46.0
x = [
    [0, 1, 0]
    [0, 1, 0]
    [1, 0, 0]
    [0, 0, 1]
    [0, 0, 0]
]

There are two sets of constraints we can dualize. It can be beneficial to apply Lagrangian Relaxation against problems composed of knapsack constraints, so we will dualize the set packing ones. The first thing we do is replace the constraints generated in lines 39-40 with a new set of variables, penalties, which take the values of the slacks on the set packing constraints. We then modify the objective function, adding Lagrangian multipliers times these penalties.

Instead of optimizing once, we do so iteratively. An important consideration is we may get nothing more than a dual bound from this process. Any integer solution is not guaranteed to be primal feasible unless it satisfies complementary slackness conditions -- for each dualized constraint either its multiplier or penalty must be zero.

Setting the intial multiplier values to 2 and using subgradient optimization with a step size of 1 / (iteration #) to adjust them, we arrive at the following model:

#!/usr/bin/env python

# This is the GAP per Wolsey, pg 182, using Lagrangian Relaxation

from gurobipy import GRB, Model

model = Model('GAP per Wolsey with Lagrangian Relaxation')
model.modelSense = GRB.MAXIMIZE
model.setParam('OutputFlag', False) # turns off solver chatter

b = [15, 15, 15]
c = [
    [ 6, 10,  1],
    [12, 12,  5],
    [15,  4,  3],
    [10,  3,  9],
    [ 8,  9,  5]
]
a = [
    [ 5,  7,  2],
    [14,  8,  7],
    [10,  6, 12],
    [ 8,  4, 15],
    [ 6, 12,  5]
]

# x[i][j] = 1 if i is assigned to j
x = []
for i in range(len(c)):
    x_i = []
    for j in c[i]:
        x_i.append(model.addVar(vtype=GRB.BINARY))
    x.append(x_i)

# As stated, the GAP has these following constraints. We dualize these into
# penalties instead, using variables so we can easily extract their values.
penalties = [model.addVar() for _ in x]
model.update()

# Dualized constraints: sum j: x_ij <= 1 for all i
for p, x_i in zip(penalties, x):
    model.addConstr(p == 1 - sum(x_i))

# sum i: a_ij * x_ij <= b[j] for all j
for j in range(len(b)):
    model.addConstr(sum(a[i][j] * x[i][j] for i in range(len(x))) <= b[j])

# u[i] = Lagrangian Multiplier for the set packing contraint i
u = [2.0] * len(x)

# Re-optimize until either we have run a certain number of iterations
# or complementary slackness conditions apply.
for k in range(1, 101):
    # max sum i,j: c_ij * x_ij
    model.setObjective(
        sum(
            # Original objective function
            sum(c_ij * x_ij for c_ij, x_ij in zip(c_i, x_i))
            for c_i, x_i in zip(c, x)
        ) + sum (
            # Penalties for dualized constraints
            u_j * p_j for u_j, p_j in zip(u, penalties)
        )
    )
    model.optimize()

    print 'iteration', k, 'obj =', model.objVal, \
        'u =', u, 'penalties =', [p.x for p in penalties]

    # Test for complementary slackness
    stop = True
    eps = 10e-6
    for u_i, p_i in zip(u, penalties):
        if abs(u_i) > eps and abs(p_i.x) > eps:
            stop = False
            break

    if stop:
        print 'primal feasible & optimal'
        break

    else:
        s = 1.0 / k
        for i in range(len(x)):
            u[i] = max(u[i] - s*(penalties[i].x), 0.0)

# Pull objective and variable values out of model
print 'objective =', model.objVal
print 'x = ['
for x_i in x:
    print '   ', [1 if x_ij.x >= 0.5 else 0 for x_ij in x_i]
print ']'

On this problem LR converges very quickly to an optimal solution, as we can see here:

iteration 1 obj = 48.0 u = [2.0, 2.0, 2.0, 2.0, 2.0] penalties = [0.0, 0.0, 0.0, 0.0, 1.0]
iteration 2 obj = 47.0 u = [2.0, 2.0, 2.0, 2.0, 1.0] penalties = [0.0, 0.0, 0.0, 0.0, 1.0]
iteration 3 obj = 46.5 u = [2.0, 2.0, 2.0, 2.0, 0.5] penalties = [0.0, 0.0, 0.0, 0.0, 1.0]
iteration 4 obj = 46.1666666667 u = [2.0, 2.0, 2.0, 2.0, 0.16666666666666669] penalties = [0.0, 0.0, 0.0, 0.0, 1.0]
iteration 5 obj = 46.0 u = [2.0, 2.0, 2.0, 2.0, 0.0] penalties = [0.0, 0.0, 0.0, 0.0, 1.0]
primal feasible & optimal
objective = 46.0
x = [
    [0, 1, 0]
    [0, 1, 0]
    [1, 0, 0]
    [0, 0, 1]
    [0, 0, 0]
]

Exercise for the reader: change the script to dualize the knapsack constraints instead of the set packing constraints. What is the result of this change in terms of convergence?

Thursday, March 29, 2012

python-zibopt 0.7.2 released

python-zibopt 0.7.2 has been released and the big news is that it supports both Python 2 and 3. (Thanks to Guillaume Sagnol for the patch!)

Other changes include:

  • Added magic square and normal magic square examples
  • Fixed issue with single-variable objectives subject to constraints:
    solver += x >= y
    solution = solver.minimize(objective=x)
  • Constant objectives no longer raise errors:
    solver.minimize(objective=3)
  • Fixed obscure constraint issue:
    solver += x <= (y-1) + 100*(1-z)

From here I'm hoping to focus my efforts on the Cython port, which should give a better basis for more advanced functionality in future versions.

Tuesday, March 27, 2012

Monday, March 26, 2012

a few thoughts on music and computer programming

Greg Linch started a discussion recently on Twitter regarding music and coding. The apex of the conversation was Zed Shaw pointing to his own well-stated comments on the matter.

I got back just over a couple hours ago from a weekly rehearsal, so I've got a number of things I'd like to mention swimming around in my head. Here are a few:

  1. To "Linguistic Skills" I'd add that musical notation is specifically built to enable certain domain-specific elements. Its staved nature allows a high degree of parallelism. I can't think of any other common language that succinctly expresses multiple concurrent paths in a readable format. One can glean a great deal of information about independently moving pieces from something like an orchestral score very quickly, in ways that we cannot with programming or other languages. Thinking concurrently is increasingly important in our contemporary computing environment.

  2. Musical notation is evolving rapidly. If one studies modern music, one will be constantly presented with new types and forms of instruction. In a way, this is like introducing new op-codes or programming constructs and can help form a desirable malleability of mind. Consider the notation, for example, in Penderecki's "Threnody for the Victims of Hiroshima" or much crazier looking things by composers such as Xenakis. Sometimes learning a single piece requires understanding the structure of an entirely new language.

  3. Music encourages team playing. In a chamber ensemble, one is expected to become so in sync with the other members that you do things like breathe together. A well-rehearsed group will practice such things as playing notes slowly to try and achieve better temperament (consistent intonation which should result in more commonly placed overtones) and vibratos that waver at the same rates. Communication is rehearsed in the form of breaths, nods, eye contact, and other physical gestures (chamber players cannot over-communicate). Most chamber music involves a lot of passing back and forth of musical phrases, which can be analogous to message passing in a synchronous computing system.

  4. Interpretation of a piece is largely subjective and therefore requires lots of compromise in order to become cohesive. This sort of compromise doesn't just mean allowing someone else to have his or her way, but means accepting and buying into it once a decision is made. If a member of my group plays a passage a certain way, I'm now expected to play that passage in a similar manner. Once a design decision is agreed upon, the impetus is on the players to accept and use it to improve the overall structure.