adventures in optimization

abusing for-loops in scala

2012-02-15T14:45:00.000-05:00

This is not usually the case, but sometimes it's nice to avoid nested iterators. If a block of code has too many, they can become visually disruptive. For an example, take a look at line 52 of this Sudoku example.

In these cases I find Python's product generator particularly useful. Consider the output of the following code:

from itertools import product
for i, j, k in product(range(2), range(10,12), range(20,22)):
    print(i, j, k)

# Output:
# 0 10 20
# 0 10 21
# 0 11 20
# 0 11 21
# 1 10 20
# 1 10 21
# 1 11 20
# 1 11 21

This behaves just like three nested for loops. You can give it an arbitrary number of sequences, too. Scala has a similar construct built into its for loops which is quite nice. Just separate your iterables by semicolons:

for (
  a <- 0 to 1;
  b <- 10 to 11;
  c <- 20 to 21
) println(a + " " + b + " " + c)

// 0 10 20
// ...

The real power here lies in the ability to apply filters to your iterables. Here we only execute the nested loops if a is even:

for (
  a <- 0 to 9 if a % 2 == 0;
  b <- 10 to 11;
  c <- 20 to 21
) println(a + " " + b + " " + c)

// 0 10 20
// 0 10 21
// 0 11 20
// 0 11 21
// 2 10 20
// 2 10 21
// 2 11 20
// 2 11 21
// 4 10 20
// ...

Of course, the more expressive a syntactic construct, the more ripe it is for abuse. It's not as blindingly obvious what the differences between the above and the following are as it would be with nesting, especially since they have the same output:

for (
  a <- 0 to 9;
  b <- 10 to 11;
  c <- 20 to 21 if a % 2 == 0
) println(a + " " + b + " " + c)

// 0 10 20
// ...

To see this a bit more clearly, we can add tests against println(...) to the innermost loop. For the uninitiated, println(...) returns Unit (think None or null), which is the same as (). The println(...) calls are essentially side effects:

for (
  a <- 0 to 9 if a % 2 == 0;
  b <- 10 to 11;
  c <- 20 to 21 if println("iterating") == ()
) println(a + " " + b + " " + c)

// iterating
// 0 10 20
// iterating
// 0 10 21
// iterating
// 0 11 20
// iterating
// 0 11 21
// iterating
// 2 10 20
// ...

We can chain our conditions nicely. Of course we have the put the println(...) first so it always executes. Note the extra iterations before a takes the value 2:

for (
  a <- 0 to 9;
  b <- 10 to 11;
  c <- 20 to 21 if println("iterating") == () if a % 2 == 0
) println(a + " " + b + " " + c)

// iterating
// 0 10 20
// iterating
// 0 10 21
// iterating
// 0 11 20
// iterating
// 0 11 21
// iterating
// iterating
// iterating
// iterating
// iterating
// 2 10 20
// ...

python-zibopt cython port started

2012-01-20T17:33:00.002-05:00

I've finally bitten the bullet and stated a branch of python-zibopt for purpose of porting it to Cython. In a sense this is like starting over, though I expect it will go fairly quickly once I get in the swing of things. The result should be a library that works in both Python 2 and 3 (currently it just supports 3) and is nearly identical to the existing API. There will probably be a few changes, but if you've only used it in the expected manner (as a modeling interface) then I doubt you'll notice any change at all.

Cython is interesting. It reminds me a bit of Inline::C, but not as terrifying. It certainly has its quirks, but its developers have done an excellent job of making the transitions between using Python and C as painless as possible. I haven't read all of the docs yet, opting for my usual method of jumping in and seeing how fast I drown, but a few points already stand out:

It isn't necessary to define every field on a C struct that you are using. You can just define the ones you use or say pass. This is extremely important, as there are so many large structs in SCIP that I thought it would be a deal-breaker.
Python strings and C char * are pretty much inter-operable. Nice.
Constants, C functions, and other things that one yanks into a .pxd (the Cython equivalent of a header file) get treated to proper Python module encapsulation and are available using dots, like foo.bar. Also nice.

So far I've just done the most basic of things: load SCIP and ask it to solve. However, there are already vast improvements in its error handling mechanisms. The C version used macros and other horrors, and it didn't catch all of SCIP's errors.

Using Cython, I was able to register an error handler with SCIP as if it were a C function (as void *, no less). I defined a simple PY_SCIP_CALL function that I wrap all my calls to SCIP in, similar to the SCIP_CALL macro that litters the SCIP source, and have it raise a Python Exception if anything is amiss. Suddenly I get all SCIP errors with appropriate text and even solver source code lines. Much better.

As a quick example, this code produces the following output. The scip.c reference in the last line shows the actual line number from ZIBopt source:

Traceback (most recent call last):
File "<string>", line 1, in <module>
File "scip.pyx", line 28, in zibopt.scip.test (src/zibopt/scip.c:606)
File "scip.pyx", line 15, in zibopt.scip.PY_SCIP_CALL (src/zibopt/scip.c:513)
Exception: [src/scip/scip.c:7500] ERROR: no node selector available

normal magic squares

2012-01-13T11:03:00.000-05:00

As a followup to yesterday's post, I created another python-zibopt example for finding Normal Magic Squares. This is similar to the Sudoku example, except that here the number of binary variables depends on the square size. In the case of Sudoku, each cell has 9 binary variables -- one for each potential value it might take. For a normal magic square, there are n^2 possible values for each cell, n^2 cells, and one variable representing the row, column, and diagonal sums. This makes a total of n^4 binary variables and one continuous variables in the model.

However, there are no big-Ms.

I think the neat part of this code is in lines 57-62. It creates sums of the n^2 variables for each cell with their appropriate coefficients (1 to n^2) and stores those expressions to make the subsequent constraint creation simpler (try doing that in your modeling language!). All made possible thanks to python-algebraic.

magic squares and big-Ms

2012-01-12T16:50:00.003-05:00

When I visited the LA PyLadies back in October of 2011, I started toying with a model for finding Magic Squares in python-zibopt. As a modeling exercise, this is fun but not too terribly challenging. Construct a square matrix of integer-valued variables and add the following constraints:

All variables >= 1.
All rows, columns, and the diagonal sum to the same value.
All variables take different values.

Admittedly, I had a few bugs to fix in the code before I could get this working. If you'd like to run it yourself, the model is here. It works against the latest development version in svn trunk of python-zibopt and python-algebraic 0.3.1. When python-zibopt 0.7.2-dev is tagged soon, it will be a part of that.

The first two constraint types are trivial to implement, and relatively easy for the solver. What I do is add a single extra variable then set it equal to the sums of each row, column, and the diagonal.

It's the third that messes things up. You can think of this as saying, for every possible pair of integer-valued variables x & y:

Either x >= y + 1 or x <= y - 1.

Why is this hard? Because we can't add both constraints to the model and maintain feasibility. What we have to do is add them in such a way that exactly one will be active for any any given solution. This requires, for each pair of variables, an additional binary variable (we will call this z) and a (possibly big) constant (M). Thus the above must be reformulated as the following before adding it to our model:

x >= (y + 1) - M*z

x <= (y - 1) + M*(1-z)

All of a sudden, here be dragons. We may not know how big or small to make M. Generally we want it as small as possible to avoid playing too much havoc with the LP relaxations of our integer programming model. It contributes to rounding errors (in the magic square problem, if I make M really big, all the variables will come back as 1). Setting M to different values may have an unpredictable effect on the solution time of a given model. So on, so forth.

Which brings us to an interesting idea:

SCIP now supports bilinear constraints out of the box. This means that I can make M a variable in the above model. (Heck, I can even make it an integer variable if I'm feeling particularly insane.)

The magic square model linked to in this post (the astute reader will notice it does not solve the normal magic square problem) provides both methods. The first command line argument it requires is the matrix size. The second one, M, is optional. If not given, it will leave M up to the solver.

I didn't expect this to perform as well as providing sensible values for M, but for small matrices it didn't perform too terribly worse either. Not quite twice the run time in most of my unscientific tests. Given the early state of MINLP development, that's pretty encouraging.

I'd love to see what one of the many far more knowledgeable OR bloggers out there has to say about this.

in defense of getting a phd

2011-12-13T11:47:00.000-05:00

Every so often I see a barrage of articles or posts discounting the value of a PhD. They make arguments along the lines of "you will spend the next six years earning minimum wage as slave to a bitter professor" or "the opportunity cost of not working until you are 30 far outweighs the benefits of such a degree".

Rarely do I see anyone argue the opposite. Granted, the fact that I am still far away from having finished a PhD at this point means that my opinion is probably worth less than the bits required to store it. However, I suspect there exists a recipe that's being overlooked by PhD naysayers.

Learn to program, learn databases, and learn the fundamentals of computer science. Learn them Really Well.
Get several years of work experience during and after college. Finish an MS or MA part time. It won't take that much longer and you'll be a much stronger candidate for jobs and PhD programs after.
Now go and do research in something else. Something that interests you. It almost doesn't matter what that is.

For me, that something else is operations research. But I've seen the formula work in completely different fields too, like political campaigns and journalism. Computing experts are scarce in most areas, and are universally necessary these days. Having multiple areas of expertise may make you the least fire-able person in your organization.

In my case, getting an MSOR qualifies me to do modeling and simulation, but if I really want to solve big, hard problems (think gas pipe lines or massive transportation networks) then I need that PhD. So much so that I sometimes consider going to school full time to speed up the process and enhance my career opportunities.

some comments on particle swarm optimization

2011-11-18T00:12:00.001-05:00

I think I'm done making changes to scala-pso for now. It finds the global minimum of Griewank's Function pretty readily, despite being somewhat magical (read: directed stochastic search). Here are a few observations, in no particular order, should anyone be interested:

PSO does a pretty decent job of finding the minimum in this case, when it has the right parameters. Make the step sizes too large, and your swarm will diverge like a herd of angry, disoriented buffalo. In this case I found that values of 0.7, 1.0, and 1.5 work well for the parameters when updating the velocity vector.
Following on the previous note, the fact that these sorts of algorithms require "metaoptimization" for their parameters is slightly distressing. Folks who are used to deterministic modeling in particular are likely to find this idea somewhat offensive.
PSO doesn't make much sense in the context of constrained optimization. To enforce boundaries, I initially just set particles on any boundary they crossed. This did not work well. Various papers suggested using the maximum distance outside of any boundary squared * 10^8 as a penalty instead. This works well, again, with the right parameters.

Here's what a swarm of 50 looks like over [-17.5, 2.5] in R^2. The global minimum is at (0, 0). I'm not including a visualization of the z-axis here -- you can find that in a previous post. Just remember that we are minimizing; think of orange as higher than blue; and larger as, well, larger. In short, we are looking for small dark blue particles. These are snapshots at the start and after successive runs of 25 swarm update iterations.

Beginning

25 iterations

50 iterations

75 iterations

100 iterations

The swarm always keeps its best found point around. In Genetic Algorithms this is known as elitism. Some algorithms use that for updating each particle's velocity. I am using the current population best, which promotes diversity a bit more. These are 95% confidence intervals of the objective (z), x1, and x2 using the same parameters, region, for 100 runs and 25 iterations (it gets kind of boring after 25):

Best objective value

Best X1

Best X2

Finally, these are 95% CIs for every particle in each of the 100 swarms, this time for 100 iterations:

Population objective value

Population X1

Population X2

griewank's function in processing

2011-11-14T00:12:00.001-05:00

I'm a little ambivalent about using Tableau for visualizing the output of this PSO project, so I turned to Processing, which is a bit nicer and fairly easy. If I can get it to create video of the particle swarm movement without too much hassle, I'll probably stick with this. Also, the image output is higher quality.

Here's what Griewank's function looks like in R^2 from [-40,10] on both axes. Processing code for doing multi-color gradients is below. The trick is to convert between Processing's idea of the origin (where 0 is at the upper left) to the human idea of that (bottom left). This is what the assignments to x1 and y1 do.

Code (sadly, no SyntaxHighlighter support yet):


// Exterior point values
float x_min = -40;
float x_max =  10;
float y_min = -40;
float y_max =  10;

// The coefficient in Griewank's function
float c = 200.0;

float griewank(float x, float y) {
  // First convert from a width x height rectangle with (0, 0)
  // at the top left to an (x_max - x_min) x (y_max - y_min) 
  // with (x_min, y_min) at the bottom left. 
  float x1 = (x/width) * (x_max-x_min) + x_min;
  float y1 = y_max - (y/height) * (y_max-y_min);
  return (x1*x1 + y1*y1) / c - (cos(x1) + cos(y1/sqrt(2))) + 1;
} 

void setup() {
  size(600, 600);
  background(0,0,0);
  noLoop();
}

void draw() {
  // Iteration 1: find the maximum.  Yes this is ugly, but it's 
  // better than storing them in an n*n array.
  float max_z = 0.0;
  for (int i = 0; i < height; i++) {
    for (int j = 0; j < width; j++) {
      float z = griewank(i, j);
      if (z > max_z)
        max_z = z;
    }
  }
  
  // Iteration 2: color in our gradient
  for (int i = 0; i < height; i++) {
    for (int j = 0; j < width; j++) {
      float z = griewank(i, j);
      color c = color(
        round(255*z/max_z), 
        round(125*z/max_z), 
        round(35*(max_z-z)/max_z)
      );      
      set(i, j, c);
    } 
  }

  save("griewank.png");
}

preliminary code for scala-pso

2011-11-10T16:42:00.001-05:00

I posted some prelimary code for PSO to github. The interface is fairly simple. Here's an example of how it works:

package scalapso.examples

import scala.math
import scalapso.Swarm

object GriewankPSO {
    val c = 200
    val lowerBounds = Array.fill(2)(-100.0)
    val upperBounds = Array.fill(2)( 100.0)
    
    // Simple implementation of Griewank's function
    def griewank(x: Array[Double]): Double = {
        var r = 0.0
        var p = 1.0
        for ((xi, i) <- x.zipWithIndex) {
            r += xi * xi / c
            p *= math.cos(xi / math.sqrt(i+1))
        }
        
        return r - p + 1
    }
    
    def main(args: Array[String]): Unit = {
        val swarm = new Swarm(griewank, lowerBounds, upperBounds, 100)
        for (i <- 0 until 10000) {
            swarm.iterate()
            println("[" + i + "] f(" + swarm.bestPoint.mkString(", ") + ") = " + swarm.bestObjective)
        }
    }
}

Lines 8 and 9 describe the feasible region. In this case we're looking in two dimensions through [-100, 100]. By changing the arguments to Array.fill (they both have to match), we could search as many dimensions as we want over any feasible region.

Lines 12 to 21 define the Griewank Function, as I posted about before. This is the most important parameter we must define for the PSO engine. As long as it has the signature Array[Double] => Double, it can be any function, which is one of the things that's attractive about this technique. (It can also be a weakness, realistically.)

In the main method, we create a swarm of 100 particles and search for 10,000 iterations. That's about it. It does reasonably well finding a point close to the global minimum of (0, 0):


[0] f(-9.907037564106162, 3.5295131972546727) = 0.845525426092787
[1] f(-9.907037564106162, 3.5295131972546727) = 0.845525426092787
[2] f(-2.75393756437488, 5.082915206876536) = 0.3345069817102082
[... snip ...]
[9998] f(-0.023389653019847856, -0.07359085931487641) = 0.0016565668346417706
[9999] f(-0.023389653019847856, -0.07359085931487641) = 0.0016565668346417706

Next step: create some fancy visualizations.

idiom vs. naivete in scala

2011-11-09T15:51:00.000-05:00

This post examines two methods for computing the Griewank Function in Scala. It compares the performance of naive versus idiomatic code for this case.

The Griewank Function is computed as follows, given an arbitrary vector of real inputs, x, and a constant, c:

$f(x) = \frac{1}{c}\sum_{i=1}^n{x_i^2} - \prod_{i=1}^n{\cos{\frac{x_i}{\sqrt{i}}}} + 1$

Obviously this will require some sort of iteration, and for my PSO library I need it extremely fast as it will be called quite a bit. A naive implementation in Scala might look something like this:

import scala.math

def griewank(x: Array[Double], c: Double): Double = {
    var r = 0.0
    var p = 1.0
    for ((xi, i) <- x.zipWithIndex) {
        r += xi * xi / c
        p *= math.cos(xi / math.sqrt(i+1))
    }
 
    return r - p + 1
}

Here we keep two accumulator variables, one for the sum and one for the product. We iterate once over x. Using x.zipWithIndex returns each value of x with its associated index, allowing us to compute the square roots in the denominator. This functions similarly to enumerate in Python, but with the order of the resulting tuples reversed.

While learning a new language, it's best to try and be as idiomatic as possible. This was a main point of the classic book, Effective Perl Programming. One can write Perl, Python, Java, etc. in any language, but when in Scala one should try and do as the Scalans do. Thus I spent a little time reformulating this implementation to use the sum and product folds that live on iterators. Here is a more Scala-esque way of doing things:

def griewank2(x: Array[Double], c: Double): Double = {
    val r2 = x.map(y => y * y / c).sum
    val p2 = x.zipWithIndex.map(
        y => math.cos(y._1 / math.sqrt(y._2 + 1))
    ).product
 
    return r2 - p2 + 1
}

For the first term, x.map takes in a lambda function and then folds it with sum. We do the same thing for the second term with zipWithIndex, using _1 and _2 to extract the tuple components and product for folding.

But there's something wrong here, which is that the latter implementation iterates twice over x. The following benchmark shows the naive implementation is about 25% faster on my machine than the idiomatic one. I guess I'll be sticking with the former.

package scalapso.examples

import scala.math
import scala.util.Random

object GriewankTest {
    // Simple implementation of Griewank's function 
    def griewank(x: Array[Double], c: Double): Double = {
        var r = 0.0
        var p = 1.0
        for ((xi, i) <- x.zipWithIndex) {
            r += xi * xi / c
            p *= math.cos(xi / math.sqrt(i+1))
        }
        
        return r - p + 1
    }
    
    // Alternative implementation using folding
    def griewank2(x: Array[Double], c: Double): Double = {
        val r2 = x.map(y => y * y / c).sum
        val p2 = x.zipWithIndex.map(
            y => math.cos(y._1 / math.sqrt(y._2 + 1))
        ).product
        
        return r2 - p2 + 1
    }
    
    def main(args: Array[String]): Unit = {
        // Generate a bunch of random points in 100 dimensions
        val rand = new Random
        val points: Array[Array[Double]] = Array.ofDim(100000, 100)
        for (point <- points)
            for (i <- 0 until point.length)
                point(i) = rand.nextDouble
    
        val c = 200
                
        // Time both implementations
        val t1 = System.currentTimeMillis
        for (point <- points)
            griewank(point, c)
        val t2 = System.currentTimeMillis    
        for (point <- points)
            griewank2(point, c)
        val t3 = System.currentTimeMillis
        
        println("simple implementation: " + (t2-t1))
        println("folded implementation: " + (t3-t2))
    }
}

Here's the output for good measure:


simple implementation: 1236
folded implementation: 1627

scala, maven, and hibernate quickstart

2011-11-07T14:06:00.000-05:00

Lately I've started a new project which will use data from pre-existing tables. It has to target the JVM, but that's its only real technical restriction. I'm taking this moment to modernize some of the components in our application stack, such as setting a better precedent for dependency handling and connecting to data tables. These are notes from the last few days of my infrastructural futzing around in Scala, Maven, and Hibernate.

In this case the tables I need to connect to already exist, so I don't need any of Hibernate's functionality for building out schema. For the sake of example, let's say a table I need to use (we'll call it Foo) looks like the following:

foo_id: integer primary key
type: varchar(5)
baz: integer
xyzzy: date

I create a new Scala project called My Awesome App using ScalaIDE and convert it to a Maven project. I then build out my directory and package structure so it resembles this:

pom.xml
src
- main.resources
  - hibernate.cfg.xml
- my.awesome.app
  - Foo.scala
  - TestFoo.scala

I update the pom.xml file, telling Maven where to get Hibernate. Normally one would put the database driver library here too, but I'm connecting to Oracle which does not, um, provide a Maven repository. So in my case I download ojbc6.jar and add it to my libraries in Eclipse, which will do for now.


    4.0.0
    my_group_id
    my-artifact-id
    0.0.1-SNAPSHOT
    My Awesome App
    
    
        
            JBoss
            JBoss Repository
            default
            http://repository.jboss.org/nexus/content/groups/public-jboss/
        
    
    
    
        
            javassist
            javassist
            3.12.1.GA
        

        
            org.hibernate
            hibernate-core
            3.6.8.Final

Now I create a Scala class representing the Foo table above. A few notes about this:

I use val instead of var, per best practices in Scala. In my case all data are read-only anyway. This makes all fields final and is a lot less annoying than using getters and setters.
Make sure to import the JPA annotations (from javax.persistence) instead of those with the same name from Hibernate. I had Eclipse import the wrong ones and figuring out what went wrong was no fun.
When a database column is named something incorrect or poor for Scala, you can give it whatever name you want on your class and use the @Column annotation to point it to the right database field.

package my.awesome.app

import java.util.Date
import javax.persistence.{Column, Entity, GeneratedValue, Id, Table}

@Entity @Table(name="Foo") class Foo {
    @Id @GeneratedValue @Column(name="foo_id") val fooId: Long = 0
    @Column(name="type") val fooType: String = ""  
    val baz: Long = 0
    val xyzzy: Date = new Date()
}

Now that I have a Scala representation of this table, I can configure Hibernate. The following is my hibernate.cfg.xml file. Note the mapping tag at the end that tells it about the annotated Foo class.





    
        
        my database driver
        my connection string
        my database user
        my database password

        
        1

        
        dialect for my database

        
        thread

        
        org.hibernate.cache.NoCacheProvider

        
        true

Finally, I create a simple main method for tesing this. Importing scala.collection.JavaConversions._ allows me to iterate over a java.util.List[Foo] instance as I would any Scala iterable.

package my.awesome.app

import org.hibernate.cfg.Configuration
import scala.collection.JavaConversions._

object TestFoo {
    def main(args: Array[String]): Unit = {
        val session = new Configuration().configure().buildSessionFactory().openSession()
        val result = session.createQuery("from Foo").list().asInstanceOf[java.util.List[Foo]]
        for (val foo <- result)
            println(foo.fooId + " | " + foo.fooType + " | " + foo.baz + " | " + foo.xyzzy)
        session.close()
    }        
}

Profit!

know your time complexities: part 2

2011-11-03T14:37:00.001-04:00

In response to this post, Ben Bitdiddle inquires:

"I understand the concept of using a companion set to remove duplicates from a list while preserving the order of its elements. But what should I do if these elements are composed of smaller pieces? For instance, say I am generating combinations of numbers in which order is unimportant. How do I make a set recognize that [1,2,3] is the same as [3,2,1] in this case?"

There are a couple points that should help here:

While lists are unhashable and therefore cannot be put into sets, tuples are perfectly capable of this. Therefore I cannot do:
```
>>> s = set()
>>> s.add([1,2,3])
Traceback (most recent call last):
 File "", line 1, in 
TypeError: unhashable type: 'list'
```
But this works just fine (extra space added for emphasis of tuple parentheses):
```
>>> s.add( (1,2,3) )
```
(3,2,1) and (1,2,3) may not hash to the same thing, but tuples are easily sortable. If I sort them before adding them to a set, they look the same:
```
>>> tuple(sorted( (3,2,1) ))
(1, 2, 3)
```

If I want to be a little fancier, I can user itertools.combinations. The following generates all unique 3-digit combinations of integers from 1 to 4:


>>> from itertools import combinations
>>> list(combinations(range(1,5), 3))
[(1, 2, 3), (1, 2, 4), (1, 3, 4), (2, 3, 4)]

Now say I want to only find those that match some condition. I can add a filter to return, say, only those 3-digit combinations of integers from 1 to 6 that multiply to a number divisible by 10:


>>> list(filter(
        lambda x: not (x[0]*x[1]*x[2]) % 10, 
        combinations(range(1, 7), 3)
    ))
[(1, 2, 5), (1, 4, 5), (1, 5, 6), (2, 3, 5), (2, 4, 5),
 (2, 5, 6), (3, 4, 5), (3, 5, 6), (4, 5, 6)]

particle swarm optimization project

2011-11-01T16:10:00.001-04:00

For my Numerical Analysis class, we've been tasked with implementing Particle Swarm Optimization (PSO) against the Griewank Function. A few notes about this:

Griewank's Function is unconstrained and non-convex. Due to the latter characteristic, its global optimum may not be easily found using traditional nonlinear optimization.
It takes in a vector of arbitrary dimension and returns a single output. The global minimum is at x = 0, f(x) = 0, making it simple to verify results. (We're looking for the global minimum, in case that wasn't clear.)
Particle Swarm "Optimization" is a bit of a misnomer. Like other types of swarm intelligence or evolutionary computation (EC), it would be better called a "directed stochastic search" than an optimization method. There's nothing wrong with that; these are really useful in problems that are non-convex or otherwise ill defined. (If this statement piques your interest, I point you to this classic paper.)

Unfortunately for me, I've about had my fill of implementing EC techniques. I actually wish I could get back some of the hours I've spent on that, but until somebody fires up a Genetic Programming engine on a Beowulf cluster to invent a time machine, that's not going to happen. Still, a few modifications should make the project really interesting:

Instead of using Matlab or a procedural language, implement it using message passing in either Erlang or Scala, where each particle has its own process or actor. This will make it scale nicely over a network.
Compare it with a simple hill climbing algorithm that starts at discretized points, or possibly with Newton's Method.
Use Tableau for visualizing the output.

Here's what Griewank's Function looks like in two dimensions. Orange corresponds to higher objective values. The global minimum is at x1 = x2 = 0, in the dark blue near the top right:

I'll post here and stick the code on github or some other place when it's coming along.

know your time complexities

2011-10-25T17:42:00.001-04:00

This is based on a lightning talk I gave at the LA PyLadies October Hackathon.

I'm actually not going to go into anything much resembling algorithmic complexity here. What I'd like to do is present a common performance anti-pattern that I see from novice programmers about once every year or so. If I can prevent one person from committing this error, this post will have achieved its goal. I'd also like to show how an intuitive understanding of time required by operations in relation to the size of data they operate on can be helpful.

Say you have a Big List of Things. It doesn't particularly matter what these things are. Often they might be objects or dictionaries of denormalized data. In this example we'll use numbers. Let's generate a list of 1 million integers, each randomly chosen from the first 100 thousand natural numbers:

import random

choices = range(100000)
x = [random.choice(choices) for i in xrange(1000000)]

Now say you want to remove (or aggregate, or structure) duplicate data while keeping them in order of appearance. Intuitively, this seems simple enough. A first solution might involve creating a new empty list, iterating over x, and only appending those items that are not already in the new list:

The Bad Way:

order = []
for i in x:
    if i not in order:
        order.append(i)

Try running this. What's wrong with it?

The issue is the conditional on line 3. In the worst case, it could look at every item in the order list for each item in x. If the list is big, as it is in our example, that wastes a lot of cycles. We can reason that we can improve the performance of our code by replacing this conditional with something faster.

Given that sets have near constant time for membership tests, one solution is to create a companion data structure, which we'll call seen. Being a set, it doesn't care about the order of the items, but it will allow us to test for membership quickly:

The Good Way:

order = []
seen = set()
for i in x:
    if i not in seen:
        seen.add(i)
        order.append(i)

Now try running this. Better?

Not that this is the best way to perform this particular action. If you aren't familiar with it, take a look at the groupby function from itertools, which is what I will sometimes reach for in a case like this.

python-algebraic

2011-08-24T00:56:00.002-04:00

After talking with a number of people and processing various bits of feedback, I've come to a few realizations regarding a shift in direction for python-zibopt.

It will be extremely painful, if not impossible, to expose solver components such as branching rules to Python using the C API. I almost got this working, but ran headlong into the problem of exception handling. Consider what happens when a branching rule or other component is written in Python and raises and exception. This exception must then be captured, handed off to the C API in python-zibopt, managed by SCIP, and then handed back to Python. Nasty stuff.
The only right way to target both Python 2 and 3 with a C extension is using Cython. This should also help with the previous item. Without porting to Cython, supporting multiple Python versions would require horrible and painful sessions of implementing unmaintainable macros.
The flexible expression handling syntax of python-zibopt may be useful to other libraries, so I am spinning it out into its own project called python-algebraic. It supports both Python 2 and 3 seamlessly. Right now it only handles those components necessary for SCIP, but with a little work it could develop into a useful library for other solver interfaces.

Lately I've felt a little like python-zibopt is taking two steps forward and one step back, but I think porting the C extension parts to Cython will help with that somewhat. It will certainly simplify the implementation.

python-zibopt 0.7 released

2011-07-07T15:52:00.005-04:00

This is mostly a documentation and bug fixing release. A few notes:

The build instructions now link python-zibopt against Ipopt for improved performance on nonlinear models.
Module docs are now being generated using reStructuredText and Sphinx. They live here.
The new expressions broke chained inequalities against a single variable in v0.6. These now work as expected: 2 <= x <= 4.
Nonlinear objective functions have to create a new variable and maximize or minimize the value of that variable. These use >= and <= in constraints for the new variable.
Dividing by variables (1/x) raises an exception.

python-zibopt module docs

2011-07-01T16:28:00.000-04:00

I needed a change of pace, so I took some time to convert all the docstrings to reStructuredText and integrate with Sphinx, the documentation builder. Apart from having to fix a few things to get Sphinx 1.1pre working in Python 3 (apparently it has to actually load a given module to generate its documentation), this was not too painful. I should have jumped on the reST bandwagon some time ago.

The module docs can be found on PyPI here.

bilinear expressions in python-zibopt 0.6

2011-06-25T12:50:00.003-04:00

The latest development release of python-zibopt supports bilinear constraints and objective functions and can be built against SCIP 2.0.1 and Python 3.2. This is very much a development version, so your mileage may vary. Algebraic expression handling has been greatly improved, along with a few other things such as the ability to remove constraints.

For 0.7 and 0.8 we are attempting to improve the build situation, providing .deb and .rpm files if possible. There will also likely be access to such properties as shadow costs and reduced costs, in addition to the ability to start writing branching rules and other solver internals directly in Python. In addition to its algebraic syntax, it is the last feature that will differentiate python-zibopt somewhat from similar libraries, in my opinion.

In 0.6 I had to completely change the way expressions are handled in order to support bilinear constraints. The following types of constraints are now possible:

solver += 2*x**2 + (3*y)**2 <= z**2 + 4
solver += (x + y) * (x - y) >= 0
solver += (x/4)**2 == 3
solver += 2 <= x <= 4

A few notes about the new features:

Expressions can exist on either side of an inequality. Previously, one side had to be a constant.
SCIP supports bilinear terms, but that is all. x**3 is invalid. Exponents must be integers from 0 to 2. Similarly, x*y is allowed but x*y*z is not.
Dividing by variables, like 1/x, is not allowed.
Chained inequalities are possible when they apply to a single variable, such as 2 <= x <= 4, but not to multiple variables as it can be impossible determine lower and upper bound constants. Consider x <= y <= z.

All constraints can be stored in Python variables using this syntax and solver.constraint(...), then removed from or added back to a solver later:

constraint = solver.constraint(2 <= x <= 4)
solver -= constraint # remove from formulation
solver += constraint # add back to formulation

When SCIP variables are subjected to mathematical operators, they are wrapped up into expression objects. These are essentially sums of terms with associated coefficients. Terms are tuples of the multiplied variables, so x becomes (x,) and x*y becomes (x, y). Thus the expression 2*x**2 + (3*y)**2 + 4*x*y - y is stored as:

{
    (x, x): 2,
    (y, y): 9,
    (x, y): 4,
    (y,): -1
}

The above requires that variables be orderable so x*y is stored the same way as y*x. Expressions handle applying algebraic logic and storing inequalities by overriding Python's mathematical functions like __add__ and __le__.

An interesting case arose when trying to deal with chained inequalities of a single variable, such as 2 <= x <= 4. Initially, I had variables return a new expression when they were subjected to <=, >= or ==. These methods normally return self, so it seemed like this might work. However, it ended up creating two expressions in this case since Python evalutes it as (2 <= x) and (x <= 4) instead of the more LISP-ish ((2 <= x) <= 4) that I was looking for. I've gotten around this in the source trunk (0.7) by having variables subclass expressions. This required adding some inelegant code to clear expression bounds off of variables once they are added to constraints, otherwise the logic for removing constraints could malfunction in single variable cases.

optimal cd laddering

2011-06-20T22:23:00.002-04:00

Here's an application that came up in conversation recently which I couldn't find any existing optimization models for: CD laddering. There are general methods for doing this, but I was interested in creating one that could take in the following inputs and produce an optimal schedule to maximize expected revenue:

A starting amount of money.
A minimum amount of money to have available each period (month) in cash and expiring CDs.
A maximum amount of money to put in any given CD in a period.
A list of CDs, their APYs, and their minimum investment amounts.

The model can be found here. It takes JSON input and makes the following basic assumptions, among other things:

APYs and investment minimums stay the same. That is, CD investment options are the same from one period to the next.
All CDs are compounded monthly.

Of course, a really useful attempt at something like this would try to predict future CD rates and simulate various scenarios to determine best and worst case. It would also deal with the fact that some CDs are compounded daily while others compounded monthly. However, given current CD rates, this probably isn't worth too much more thought at this point.

python-zibopt google group

2011-06-20T12:42:00.004-04:00

There is a new Google group intended for discussion of python-zibopt announcements, usage, and features. Please consider joining if you use this API at all.

syntax highlighting for r

2011-06-13T11:44:00.002-04:00

Thanks to Paul Rubin and his blog for pointing me to Pretty R. I've been using SyntaxHighlighter for displaying Python and C, but it doesn't have a brush for highlighting R code. Pretty R is a nice substitute for the interim, as you can see in this updated post.

a few notes on deterministic vs. stochastic simulation

2011-06-09T15:49:00.002-04:00

I find I have to build simulations with increasing frequency in my work and life. Usually this indicates I'm faced with one of the following situations:

The need for a quick estimate regarding the quantitative behavior of some situation.
The desire to verify the result of a computation or assumption.
A situation which is too complex or random to effectively model or understand.

Anyone familiar at all with simulation will recognize the last item as the motivating force of the entire field. Simulation models tend to take over when systems become so complex that understanding them is prohibitive in cost and time or entirely infeasible. In a simulation, the modeler can focus on individual interactions between entities while still hoping for useful output in the form of descriptive statistics.

As such, simulations are nearly always stochastic. The output of a simulation, whether it be the mean time to service upon entering a queue or the number of fish alive in a pond, is determined by a number of random inputs. It is estimated by looking at a sample of the entire, often infinite, problem space and therefore must be described in terms of mean and variance.

For me, simulation building usually follows a process roughly like this:

Work with a domain expert to understand the process under study.
Convert this process into a deterministic simulation (no randomness).
Verify the output of the deterministic simulation.
Anlyze the inputs of the simulation to determine their probability distributions.
Convert the deterministic simulation to a stochastic simulation.

The reason for creating a simulation without randomness first is that it can be difficult or impossible to verify its correctness otherwise. Thus one may focus on the simulation logic first before analyzing and adding sources of randomness.

Where the procedure breaks down is after the third step. Domain experts are often happy to share their knowledge about systems to aid in designing simulations, and typically can understand the resulting abstractions. They are also invaluable in verifying simulation output. However, they are unlikely to understand why it is necessary to add randomness to a system that they already perceive as functional. Further, doing so can be just as difficult and time consuming as the initial model development and therefore requires justification.

This can be a quandary for the model builder. How does one communicate the need to incorporate randomness to decision makers who lack understanding of probability? It is trivially easy to construct simulations that use the same input parameters but yield drastically different outputs. Consider the code below, which simulates two events occurring and counts the number of times event b happens before event a.

import random

def sim_stochastic(event_a_lambda, event_b_lambda):
    # Returns 0 if event A arrives first, 1 if event B arrives first

    # Calculate next arrival time for each event randomly.
    event_a_arrival = random.expovariate(event_a_lambda)
    event_b_arrival = random.expovariate(event_b_lambda)

    return 0.0 if event_a_arrival <= event_b_arrival else 1.0

def sim_deterministic(event_a_lambda, event_b_lambda):
    # Returns 0 if event A arrives first, 1 if event B arrives first

    # Calculate next arrival time for each event deterministically.
    event_a_arrival = 1.0 / event_a_lambda
    event_b_arrival = 1.0 / event_b_lambda

    return 0.0 if event_a_arrival <= event_b_arrival else 1.0

if __name__ == '__main__':
    event_a_lambda = 0.3
    event_b_lambda = 0.5

    repetitions = 10000

    for sim in (sim_stochastic, sim_deterministic):
        output = [
            sim(event_a_lambda, event_b_lambda)
            for _ in range(repetitions)
        ]
        event_b_first = 100.0 * (sum(output) / len(output))
        print('event b is first %0.1f%% of the time' % event_b_first)

Both simulations use the same input parameter, but the second one is essentially wrong as b will always happen first. In the stochastic version, we use exponential distributions for the inputs and obtain an output that verifies our basic understanding of these distributions.

event b is first 63.0% of the time
event b is first 100.0% of the time

How about you? How do you discuss the need to model a random world with decision makers?

new in python-zibopt: constraint removal

2011-05-24T17:47:00.002-04:00

At long last, python-zibopt has constraint removal and a nice interface to go with it. The idea here is that you store constraints in variables once they are created and can remove them from the solver, and add them back in, as you see fit.

The feature is in trunk, so it will be part of 0.6 when that is released. A few notes about using it:

Constraint instances are only returned using solver.constraint(...). The algebraic syntax can't actually return the constraints, so that limits you from removing them later.
It is possible to reinstall a removed constraint using +=. See below.
Removal or addition of constraints necessitates restarting the solver.

Here's an example. I have to admit, I'm getting pretty excited about this release.

from zibopt import scip

solver = scip.solver()

x1 = solver.variable(scip.INTEGER)
x2 = solver.variable(scip.INTEGER)

c1 = solver.constraint(upper=4, coefficients={x1:2, x2:2})
c2 = solver.constraint(upper=3, coefficients={x1:2, x2:2})

print(solver.maximize(objective=x1+x2).objective) # 1.0

solver -= c2
print(solver.maximize(objective=x1+x2).objective) # 2.0

solver += c2
print(solver.maximize(objective=x1+x2).objective) # 1.0

python-zibopt works with python 3.2 in trunk

2011-05-21T16:05:00.005-04:00

python-zibopt is now ported to Python 3.2 in the repository trunk. All the tests pass and the examples are working, so if anybody wants to get latest and play around, that should be OK. Some notes about the process of porting the C extension, in case anyone else finds them useful:

Changing special method names in a duck typed language makes porting painful. For instance, __nonzero__ in Python 2 became
__bool__ in Python 3. Without proper unit tests, I may never even have caught this, since Python will go ahead and provide default implementations for these sorts of methods, completely ignoring the existence of older versions of them.

The easiest part was getting the examples to work. These just required changing of print statements and adding calls to list around functions that have become lazy (zip, for instance). This is a good sign for folks that have existing model code in python-zibopt.

init_foo should now by PyInit_foo. Further, each C type needs a PyModuleDef and should use PyModule_Create instead of Py_InitModule3. The current docs show this accurately.

It's difficult to find examples of certain, rather fundamental things in the documentation. Particularly, I had to pull the right method for overriding __getattr__ and __setattr__ in C from the code from pyexpat. Since strings in Python are all now unicode, what once looked like:

static PyObject* branching_rule_getattr(branching_rule *self, PyObject *attr_name) {
   char *attr;

   // Check and make sure we have a string as attribute name...
   if (PyString_Check(attr_name)) {
       attr = PyString_AsString(attr_name);

       if (!strcmp(attr, "maxbounddist"))
           return Py_BuildValue("d", SCIPbranchruleGetMaxbounddist(self->branch));
       if (!strcmp(attr, "maxdepth"))
           return Py_BuildValue("i", SCIPbranchruleGetMaxdepth(self->branch));
       if (!strcmp(attr, "priority"))
           return Py_BuildValue("i", SCIPbranchruleGetPriority(self->branch));
   }
   return PyObject_GenericGetAttr(self, attr_name);
}

now requires use of PyUnicode_Check and PyUnicode_CompareWithASCIIString (note the == 0 after these function calls):

static PyObject* branching_rule_getattr(branching_rule *self, PyObject *attr_name) {
   // Check and make sure we have a string as attribute name...
   if (PyUnicode_Check(attr_name)) {
       if (PyUnicode_CompareWithASCIIString(attr_name, "maxbounddist") == 0)
           return Py_BuildValue("d", SCIPbranchruleGetMaxbounddist(self->branch));
       if (PyUnicode_CompareWithASCIIString(attr_name, "maxdepth") == 0)
           return Py_BuildValue("i", SCIPbranchruleGetMaxdepth(self->branch));
       if (PyUnicode_CompareWithASCIIString(attr_name, "priority") == 0)
           return Py_BuildValue("i", SCIPbranchruleGetPriority(self->branch));
   }
   return PyObject_GenericGetAttr((PyObject *) self, attr_name);
}

Changes to the __setattr__ code are similar.

SCIP provides an opt-shared makefile for Linux x86_64, but SoPlex does not, for whatever reason. Here is a diff to create a make.linux.x86_64.gnu.opt-shared file from make.linux.x86_64.gnu.opt in SoPlex 1.5.0 in case you need that. The existing build instructions are generally accurate for the new version on SCIP and Python 3.2, though an update will arrive with v0.6. Your mileage may vary for now.

4,5c4,5
< CXXFLAGS= -O3 -m64 -mtune=native -fomit-frame-pointer # -malign-double -ffast-math < LDFLAGS = -lm -static --- > CXXFLAGS= -O3 -m64 -mtune=native -fomit-frame-pointer -fPIC # -malign-double -ffast-math
> LDFLAGS = -lm -Wl,-rpath,$(CURDIR)/$(LIBDIR)
8a9,15
> LIBBUILD        =       $(CXX)
> LIBBUILDFLAGS   = -shared -FPIC
> LIBBUILD_o      =       -o # the trailing space is important
> LIBEXT  = so
> ARFLAGS         =
> RANLIB          =
>

joy in the time of the python futures module

2011-05-19T07:42:00.000-04:00

Starting with a does of realism, it's possible this will turn out like the day when coroutines were introduced into Python 2.5. At the time I was extremely excited. Upon hearing the news I spent several hours trying to convince developers at my then-employer that the only way we could survive as an organization was to immediately abandon all our existing Java infrastructure, porting it to a new and beautiful world based on finite state machines implemented using Python coroutines. After a day of hand waving over a proof of concept, we all continued about our lives. Soon after, I left for a Python shop, but in the next half decade I still never found a good place to use this beloved feature in solving the daily challenges of my professional life.

As I come to terms more with switching to Python 3.2, the futures module emerges as a source of similar excitement. This is one of those would-have-made-my-life-so-much-easier features I wish I'd had years ago, and is almost reason in itself to upgrade from Python 2.7. Who cares if none of your libraries have been ported yet?

I think the real strength of this library springs from its ability to take any pre-existing function and distribute it over a process pool. Here is an example that computes minimum spanning trees for fully connected graphs. For purposes of testing, we generate 100-node fully connected graphs with random arc weights between 1 and 10. We then find their minimum spanning trees using the pygraph library, using groups of 4, 8, ..., 28 random graphs in serial and in parallel. Note how easy it is to take the minimum spanning tree function and map it over a process pool without any changes to its code.

from concurrent import futures
from csv import writer
from pygraph.algorithms import minmax
from pygraph.classes.graph import graph
import random
import string
import sys
import time

def generate_graph():
    # Generates a randomly weighted, fully-connected, undirected graph
    nodes = [str(i) for i in range(100)]

    g = graph()
    for n in nodes:
        g.add_node(n)

    # Build an edge from each node to every other node
    for i, n in enumerate(nodes):
        for o in nodes[i+1:]:
            weight = random.uniform(1, 10)
            g.add_edge((n, o), weight)

    return g

def serial_test(graphs):
    for g in graphs:
        tree = minmax.minimal_spanning_tree(g)

def parallel_test(graphs, max_workers):
    with futures.ProcessPoolExecutor(max_workers=max_workers) as executor:
        for tree in executor.map(minmax.minimal_spanning_tree, graphs):
            pass # normally we'd do something with this...

if __name__ == '__main__':
    out = writer(sys.stdout)
    out.writerow(['num graphs', 'serial time', 'parallel time'])

    # Run with a number of different randomly generated graphs
    for num_graphs in (4, 8, 12, 16, 20, 24, 28):
        graphs = [generate_graph() for _ in range(num_graphs)]

        start = time.clock()
        serial_test(graphs)
        serial_time = time.clock() - start

        start = time.clock()
        parallel_test(graphs, 4)
        parallel_time = time.clock() - start

        out.writerow([num_graphs, serial_time, parallel_time])

The output of this script shows that we get a fairly linear speedup in this particular example with little to no effort.

Given that the box I'm running this on has 4 cores, it's a little odd that the speedup factor is more like 2. It's probably just that the machine has a lot going on, so it's not really worth investigating right now. At the very least, each core is kept busy when the test forks.

python-zibopt upgrade on the way

2011-05-18T16:18:00.002-04:00

To anyone out there who is using python-zibopt, I have begun the painful processing of upgrading it to both Python 3.2 and SCIP 2.0.1 simultaneously. I think it makes sense to do this now, because both C APIs have changed and I'd rather just have to deal with this once.

There is also a compelling reason to upgrade the Python API now in PEP 384: Defining a Stable ABI. Hopefully, the Python C API should remain stable after this and thus not require future changes for C extension maintainers. It's hard enough keeping up with changes to one major system, let alone the underlying interpreter.

With any luck I'll get this working in the next few weeks, along with the new nonlinear constraints available in SCIP. That will be packaged into a 0.6 release, then I can start to focus on some of the remaining issues such as figuring out a better method for building the library.

rip bernard greenhouse

2011-05-16T19:37:00.003-04:00

R.I.P. Bernard Greenhouse. As with so many other cellists of my generation, you were my teacher's teacher. While I never saw you perform in person, my lessons were often sprinkled with references to you. You will be missed, and it is lost on no one that you passed three weeks after you stopped practicing.

thoughts on mandarin in the media

2011-05-16T16:15:00.002-04:00

Thanks, MacDiva, for posting a link to this article.

I believe this is a good example of the recent media movement trying to shock Americans into believing we'll fall behind as a culture if we don't all go out and learn Mandarin. While it's very important to include second (or third) languages in the education of our children (and I desparately wish I'd taken something other than Latin in school), there's a lot that's disingenuous about the sentiments expressed here. A few points:

If person A learns the language of person B, it is not required for person B to reciprocate in order to establish effective communication.
The impact of having every young Chinese citizen learn English will hopefully be that most of them can speak English.
If this program is successful, there will come a time when native English speakers fluent in Mandarin are no longer necessary in such industries as tourism, trade, and international relations in China. Those functions will be more often assigned to Chinese citizens. Thus, learning Mandarin at this time has limited utility.
When a large body of people devotes effort to learning a foreign language, that is a recognition of the importance of that language. If anything, this should assist in the established dominance of English as an international medium for communication.

Perhaps I'm wrong in my thinking, but I don't see this as a good reason to develop a national, or even individual, effort to learn Mandarin. One could do so entirely reasonably for personal betterment and cultural exposure, but this infographic looks like media scaremongering to me.

zen of lisp / lisp koans

2011-04-28T18:41:00.008-04:00

There used to be a website of collected mythology about LISP masters and students, mostly at MIT during the 60s and 70s. I believe it was called the Zen of LISP or something to that effect, but can no longer seem to find it. After an extensive search, I've resigned myself to pulling as many of these as I can find out of my local fortune file. Here are a few gems:

A novice was trying to fix a broken lisp machine by turning the power off and on. Knight, seeing what the student was doing spoke sternly, "You cannot fix a machine by just power-cycling it with no understanding of what is going wrong." Knight turned the machine off and on. The machine worked.

Of course, it's entirely possible to fix a machine by power-cycling if you do understand what is going on.

A famous Lisp Hacker noticed an Undergraduate sitting in front of a Xerox 1108, trying to edit a complex Klone network via a browser. Wanting to help, the Hacker clicked one of the nodes in the network with the mouse, and asked "what do you see?" Very earnestly, the Undergraduate replied, "I see a cursor." The Hacker then quickly pressed the boot toggle at the back of the keyboard, while simultaneously hitting the Undergraduate over the head with a thick Interlisp Manual. The Undergraduate was then Enlightened.

Somehow I doubt it was ever that easy. When I hit myself over the head with CLTL2E, it just hurts.

A student, in hopes of understanding the Lambda-nature, came to Greenblatt. As they spoke a Multics system hacker walked by. "Is it true", asked the student, "that PL-1 has many of the same data types as Lisp?" Almost before the student had finished his question, Greenblatt shouted, "FOO!", and hit the student with a stick."

Not entirely sure the right way to approach this one. Does it have to do with dynamic typing, metasyntactic variables, both, something else entirely?

There's one in particular I have in mind that has to do with garbage collection and cycles in reference counting, but I haven't found it yet. If you happen to know of it, please send it along.

affine scaling in r

2011-04-27T07:30:00.002-04:00

I recently stumbled across an implementation of the affine scaling interior point method for solving linear programs that I'd coded up in R once upon a time. I'm posting it here in case anyone else finds it useful. There's not a whole lot of thought given to efficiency or numerical stability, just a demonstration of the basic algorithm. Still, sometimes that's exactly what one wants.

solve.affine <- function(A, rc, x, tolerance=10^-7, R=0.999) {
  # Affine scaling method
  while (T) {
    X_diag <- diag(x)
 
    # Compute (A * X_diag^2 * A^t)-1 using Cholesky factorization.
    # This is responsible for scaling the original problem matrix.
    q <- A %*% X_diag**2 %*% t(A)
    q_inv <- chol2inv(chol(q))
 
    # lambda = q * A * X_diag^2 * c
    lambda <- q_inv %*% A %*% X_diag^2 %*% rc
 
    # c - A^t * lambda is used repeatedly
    foo <- rc - t(A) %*% lambda
 
    # We converge as s goes to zero
    s <- sqrt(sum((X_diag %*% foo)^2))
 
    # Compute new x
    x <- (x + R * X_diag^2 %*% foo / s)[,]
 
    # If s is within our tolerance, stop.
    if (abs(s) < tolerance) break
  }
  x
}

This function accepts a matrix A which contains all technological coefficients for an LP, a vector rc containing its reduced costs, and an initial point x interior to the LP's feasible region. Optional arguments to the function include a tolerance, for detecting when the method is within an acceptable distance from the optimal point, and a value for R, which must be strictly between 0 and 1 and controls scaling.

The method works by rescaling the matrix A around the current solution x. It then computes a new x such that it remains feasible and interior, which is why R cannot be 0 or 1. It requires a feasible interior point to start and only projects to other feasible interior points, so the right hand side of the LP is not required (it is implicit from the starting point). The shadow prices for each iteration are captured in the vector lambda, so the gap between primal and dual solutions is easy to compute.

We run this function against a 3x3 LP with a known solution:

max z = 5x1 + 4x2 + 3x3
st      2x1 + 3x2 +  x3 <=  5 
        4x1 +  x2 + 2x3 <= 11 
        3x1 + 4x2 + 2x3 <=  8 
        x1, x2, x3 >= 0

The optimal solution to this LP is:

z  = 13 
x1 =  2 
x2 =  0 
x3 =  1

This problem can be run against the affine scaling function by defining A with all necessary slack variables, and using an arbitrary feasible interior point:

A  <- matrix(c(
  2,3,1,1,0,0, 
  4,1,2,0,1,0, 
  3,4,2,0,0,1
), nrow=3, byrow=T) 
rc <- c(5, 4, 3, 0, 0, 0) 
x  <- c(0.5, 0.5, 0.5, 2, 7.5, 3.5) 

solution <- solve.affine(A, rc, x) 
print(solution) 
print(sum(solution * rc))

This provides an output vector that is very close to the optimal primal solution shown above. Since interior point methods converge asymptotically to optimal solutions, it is important to note that we can only ever get (extremely) close to our final optimal objective and decision variable values.

> print(solution) 
[1] 1.999998e+00 4.268595e-07 1.000002e+00 1.280579e-06 1.000005e+00 
[6] 1.280579e-06 
 
> print(sum(solution * rc))
[1] 13.00000

reformed japhs in python: scheme to python compilation

2011-04-20T11:00:00.000-04:00

I believe this is the final JAPH in this series. I actually didn't have the heart to obfuscate it. It starts with a Scheme program that prints 'just another scheme hacker', tokenizes it, parses the token stream, compiles that into Python 3.2, and executes the resulting string. If anybody else wants to obfuscate it, be my guest. Otherwise I'll just let it speak for itself.

import re

def tokenize(input):
    '''Tokenizes an input stream into a list of recognizable tokens'''
    token_res = (
        r'\(',      # open paren -> starts expression
        r'\)',      # close paren -> ends expression
        r'"[^"]*"', # quoted string (don't support \" yet)
        r'[\w?]+'   # atom
    )
    return re.findall(r'(' + '|'.join(token_res) + ')', input)

def parse(stream):
    '''Parses a token stream into a syntax tree'''
    if not stream:
        return []

    else:
        # Build a list of arguments (possibly expressions) at this level
        args = []
        while True:
            # Get the next token
            try:
                x = stream.pop(0)
            except IndexError:
                return args

            # ( and ) control the level of the tree we're at
            if x == '(':
                args.append(parse(stream))
            elif x == ')':
                return args
            else:
                args.append(x)

def compile(tree):
    '''Compiles an Scheme Abstract Syntax Tree into near-Python'''
    def compile_expr(indent, expr):
        indent += 1

        lines = [] # these will have [(indent, statement), ...] structure
        while expr:
            # Two options: expr is a string like "'" or it is a list
            if isinstance(expr, str):
                return [(
                    indent,
                    expr.replace('scheme', 'python').replace('\n', '\\n')
                )]

            else:
                start = expr.pop(0)

                if start == 'define':
                    signature = expr.pop(0)
                    lines.append((indent,
                        'def %s(%s):' % (
                            signature[0],
                            ', '.join(signature[1:])
                        )
                    ))
                    while expr:
                        lines.extend(compile_expr(indent, expr.pop(0)))

                elif start == 'if':
                    # We don't support multi-clause conditionals yet
                    clause = compile_expr(indent, expr.pop(0))[0][1]
                    lines.append((indent, 'if %s:' % clause))

                    if_true_lines = compile_expr(indent, expr.pop(0))
                    if_false_lines = compile_expr(indent, expr.pop(0))

                    lines.extend(if_true_lines)
                    lines.append((indent, 'else:'))
                    lines.extend(if_false_lines)

                elif start == 'null?':
                    # Only supports conditionals of the form (null? foo)
                    if isinstance(expr[0], str):
                        condition = expr.pop(0)
                    else:
                        condition = compile_expr(indent, expr.pop(0))[0][1]
                    return [(indent, 'not %s' % condition)]

                elif start == 'begin':
                    # This is just a series of statements, so don't indent
                    while expr:
                        lines.extend(compile_expr(indent-1, expr.pop(0)))

                elif start == 'display':
                    arguments = []
                    while expr:
                        arguments.append(
                            compile_expr(indent, expr.pop(0))[0][1]
                        )
                    lines.append((
                        indent,
                        "print(%s, end='')" % (', '.join(arguments))
                    ))

                elif start == 'car':
                    lines.append((indent, '%s[0]' % expr.pop(0)))

                elif start == 'cdr':
                    lines.append((indent, '%s[1:]' % expr.pop(0)))

                elif start == 'list':
                    arguments = []
                    while expr:
                        arguments.append(
                            compile_expr(indent, expr.pop(0))[0][1]
                        )
                    lines.append((indent, '[%s]' % ', '.join(arguments)))

                else:
                    # Assume this is a function call
                    arguments = []
                    while expr:
                        arguments.append(
                            compile_expr(indent, expr.pop(0))[0][1]
                        )
                    lines.append((
                        indent,
                        "%s(%s)" % (start, ', '.join(arguments))
                    ))

        return lines

    return [compile_expr(-1, expr) for expr in tree]

if __name__ == '__main__':
    scheme = '''
        (define (output x)
            (if (null? x)
                ""
                (begin (display (car x))
                       (if (null? (cdr x))
                           (display "\n")
                           (begin (display " ")
                                  (output (cdr x)))))))
        (output (list "just" "another" "scheme" "hacker"))
    '''
    python = ''
    for expr in compile(parse(tokenize(scheme))):
        python += '\n'.join([(' ' * 4 * x[0]) + x[1] for x in expr]) + '\n\n'
    exec(python)

reformed japhs in python: turing machine

2011-04-18T08:00:00.000-04:00

This JAPH constructs a Turing machine in order to achieve its goal. The machine accepts any string that ends in '\n' and, to assist in our purposes, allows side effects. This lets us print the value of the tape as it encounters each character. While the idea of using lambda functions as side effects in a Turing machine is a little bizarre on many levels, we work with what we have. And Python is multi-paradigmatic, so what the heck.

import re

def turing(tape, transitions):
    # The tape input comes in as a string.  We approximate an infinite
    # length tape via a hash, so we need to convert this to {index: value}
    tape_hash = {i: x for i, x in enumerate(tape)}

    # Start at 0 using our transition matrix
    index = 0
    state = 0
    while True:
        value = tape_hash.get(index, '')

        # This is a modified Turing machine: it uses regexen
        # and has side effects.  Oh well, I needed IO.
        for rule in transitions[state]:
            regex, next, direction, new_value, side_effect = rule
            if re.match(regex, value):
                # Terminal states
                if new_value in ('YES', 'NO'):
                    return new_value

                tape_hash[index] = new_value
                side_effect(value)
                index += direction
                state = next
                break

assert 'YES' == turing('just another python hacker\n', [
    # This Turing machine recognizes the language of strings that end in \n.

    # Regex rule, next state, left/right = -1/+1, new value, side effect.
    [ # State 0:
        [r'^[a-z ]$', 0, +1, '', lambda x: print(x, end='')],
        [r'^\n$', 1, +1, '', lambda x: print(x, end='')],
        [r'^.*$', 0, +1, 'NO', None],
    ],
    [ # State 1:
        [r'^$', 1, -1, 'YES', None]
    ]
])

Obfuscation again consists of converting the above code into lambda functions using Y combinators. This is a pretty fantastic programming exercise, so I've left it out of this post in case anyone wants to try. And of course the Turing machine has to return 'YES' to indicate that it accepts the string, thus the assertion. Our final obfuscated JAPH looks like this, amazingly in a single expression:

assert'''YES'''==(lambda g:(lambda f:g(lambda arg:f(f)(arg)))(lambda f:g(
lambda arg: f(f)(arg))))(lambda f: lambda q:[(lambda g:(lambda f:g(lambda
arg:f(f)(arg)))(lambda f: g(lambda arg:f(f)(arg))))(lambda f: lambda x:(x
[0][0]if x[0] and __import__('re').match(x[0][0][0],x[1])else f([x[0][1:]
,x[1]]))) ([q[3][q[1]],q[2].get(q[0],'')])[4](q[2].get(q[0],'')), (lambda
g:(lambda f:g(lambda arg:f(f)(arg))) (lambda f:g(lambda arg:f(f)(arg))))(
lambda f:lambda x:(x[0][0]if x[0] and __import__('re').match(x[0][0][0],x
[1])else f([x[0][1:],x[1]])))([q[3][q[1]],q[2].get(q[0],'')])[3]if(lambda
g:(lambda f:g(lambda arg:f(f)(arg))) (lambda f:g(lambda arg:f(f)(arg))))(
lambda f:lambda x:(x[0][0]if x[0]and __import__('re').match(x[0][0][0],x[
1]) else f([x[0][1:],x[1]])))([q[3][q[1]],q[2].get(q[0],'')])[3]in('YES',
'NO')else f([q[0]+(lambda g:(lambda f:g(lambda arg:f(f)(arg)))(lambda f:g
(lambda arg:f(f)(arg))))(lambda f:lambda x:(x[0][0]if x[0]and __import__(
're').match(x[0][0][0],x[1])else f([x[0][1:], x[1]])))([q[3][q[1]], q[2].
get(q[0],'')])[2],(lambda g:(lambda f:g(lambda arg: f(f)(arg)))(lambda f:
g(lambda arg:f(f)(arg))))(lambda f:lambda x:(x[0][0]if x[0]and __import__
('re').match(x[0][0][0],x[1])else f([x[0][1:], x[1]])))([q[3][q[1]],q[2].
get(q[0],'')])[1],q[2],q[3]])][1])([0,0,{i:x for i,x in enumerate('just '
'another python hacker\n')}, [[[r'^[a-z ]$',0,+1,'',lambda x:print(x,end=
'')], [r'^\n$',1,+1,'',lambda x:print(x, end='')],[r'^.*$',0,+1,'''NO''',
lambda x:None]], [[r'''^$''',+1,-1,'''YES''', lambda x: None or None]]]])

reformed japhs in python: huffman coding

2011-04-14T10:30:00.000-04:00

At this point, tricking python into printing strings via ever more pernicious mechanisms got a little boring. So I switched to obfuscating fundamental computer science algorithms. Here's a JAPH that takes in a Huffman coded version of 'just another python hacker', decodes, and prints it.

# Build coding tree
def build_tree(scheme):
    if scheme.startswith('*'):
        left, scheme = build_tree(scheme[1:])
        right, scheme = build_tree(scheme)
        return (left, right), scheme
    else:
        return scheme[0], scheme[1:]

def decode(tree, encoded):
    ret = ''
    node = tree
    for direction in encoded:
        if direction == '0':
            node = node[0]
        else:
            node = node[1]
        if isinstance(node, str):
            ret += node
            node = tree
    return ret

tree = build_tree('*****ju*sp*er***yct* h**ka*no')[0]
print(
    decode(tree, bin(10627344201836243859174935587).lstrip('0b').zfill(103))
)

The decoding tree is built like a true LISP-style tree as a sequence of pairs. '*' represents a branch in the tree while other characters are leaf nodes. This looks like ((((('j', 'u'), ('s', 'p')), ('e', 'r')), ((('y', 'c'), 't'), (' ', 'h'))), (('k', 'a'), ('n', 'o'))) after it's constructed.

The actual Huffman coded version of our favorite string looks like 0000000001000100101011010111011101010111001000110110000110100001010111111110011001111010100110000100011, which in base-2 encoding represents around a 50% savings in space. Well worth all the effort, right?

There's a catch here, which is that this is hard to obfuscate unless we turn it into a single expression. This means that we have to convert build_tree and decode into lambda functions. Unfortunately, they are recursive and lambda functions don't do that easily. Fortunately, we can use Y combinators and then the rest is simple. These are worth some study since they will pop up again later.

Y = lambda g: (
    lambda f: g(lambda arg: f(f)(arg))) (lambda f: g(lambda arg: f(f)(arg))
)

build_tree = Y(
    lambda f: lambda scheme: (
        (f(scheme[1:])[0], f(f(scheme[1:])[1])[0]),
        f(f(scheme[1:])[1])[1 ]
    ) if scheme.startswith('*') else (scheme[0], scheme[1:])
)

decode = Y(lambda f: lambda x: x[3]+x[1] if not x[2] else (
    f([x[0], x[0], x[2], x[3]+x[1]]) if isinstance(x[1], str) else (
        f([x[0], x[1][0], x[2][1:], x[3]]) if x[2][0] == '0' else (
            f([x[0], x[1][1], x[2][1:], x[3]])
        )
    )
))

tree = build_tree('*****ju*sp*er***yct* h**ka*no')[0]
print(
    decode([
        tree,
        tree,
        bin(10627344201836243859174935587).lstrip('0b').zfill(103), ''
    ])
)

The final version is really just a condensed (and expanded, weirdly) version of the above (again, make sure to use Python 3.2):

print((lambda t,e,s:(lambda g:(lambda f:g(lambda arg:f(f)(arg)))(lambda f:
g(lambda arg: f(f)(arg))))(lambda f:lambda x: x[3]+x[1]if not x[2]else f([
x[0],x[0],x[2],x[3]+x[1]])if isinstance(x[1],str)else f([x[0],x[1][0],x[2]
[1:],x[3]])if x[2][0]=='0'else f([x[0],x[1][1],x[2][1:],x[3]]))([t,t,e,s])
)((lambda g:(lambda f:g(lambda arg:f(f)(arg)))(lambda f:g(lambda arg:f(f)(
arg))))(lambda f:lambda p:((f(p[1:])[0],f(f(p[1:])[1])[0]),f(f(p[1:])[1])[
1])if p.startswith('*')else(p[0],p[1:]))('*****ju*sp*er***yct* h**ka*no')[
0],bin(10627344201836243859179756385-4820798).lstrip('0b').zfill(103),''))

reformed japhs in python: rolling effect

2011-04-11T09:20:00.000-04:00

Here's a JAPH composed solely for effect. For each letter in 'just another python hacker' it loops over each the characters ' abcdefghijklmnopqrstuvwxyz', printing each. Between characters it pauses for 0.05 seconds, backing up and moving on to the next if it hasn't reached the desired one yet. This achieves a sort of rolling effect by which the final string laboriously appears on our screen.

import string
import sys
import time

letters = ' ' + string.ascii_lowercase
for l in 'just another python hacker':
    for x in letters:
        print(x, end='')
        sys.stdout.flush()
        time.sleep(0.05)

        if x == l:
            break
        else:
            print('\b', end='')

print()

Locating and printing of each letter in the string is done via a list comprehension. At the end we actually have an extra line of code (the eval statement) that gives us our newline.

[[(lambda x,l:str(print(x,end=''))+str(__import__(print.
__doc__[print.__doc__.index('stdout') - 4:print.__doc__.
index('stdout')-1]).stdout.flush()) + str(__import__(''.
join(reversed('emit'))).sleep(0o5*1.01/0x64))+str(print(
'\b',end='\x09'.strip())if x!=l else'*&#'))(x1,l1)for x1
in('\x20'+getattr(__import__(type('phear').__name__+'in'
'g'),dir(__import__(type('snarf').__name__+'ing'))[13]))
[:('\x20'+getattr(__import__(type('smear').__name__+'in'
'g'),dir(__import__(type('slurp').__name__+'ing'))[13]))
.index(l1)+1]]for l1 in'''just another python hacker''']
eval('''\x20\x09eval("\x20\x09eval('\x20 print()')")''')

reformed japhs in python: rot13

2011-04-06T07:30:00.002-04:00

No series of JAPHs would be complete without ROT13. This is the example through which aspiring Perl programmers learn to use tr and its synonym y. In Perl the basic ROT13 JAPH starts as:

$foo = 'whfg nabgure crey unpxre';
$foo =~ y/a-z/n-za-m/;
print $foo;

Python has nothing quite so elegant in its default namespace. However, this does give us the opportunity to explore a little used aspect of strings: the translate method. If we construct a dictionary of ordinals we can accomplish the same thing with a touch more effort.

import string

table = {
    ord(x): ord(y) for x, y in zip(
        string.ascii_lowercase,
        string.ascii_lowercase[13:] + string.ascii_lowercase
    )
}

print('whfg nabgure clguba unpxre'.translate(table))

We obfuscate the construction of this translation dictionary and, for added measure, use getattr to find the print function off of __builtins__. This will likely only work in Python 3.2, since the order of attributes on __builtins__ matters.

getattr(vars()[list(filter(lambda _:'\x5f\x62'in _,dir
()))[0]], dir(vars()[list(filter(lambda _:'\x5f\x62'in
_, dir()))[0]])[list(filter(lambda _:_ [1].startswith(
'\x70\x72'),enumerate(dir(vars()[list(filter(lambda _:
'\x5f\x62'in _,dir()))[0]]))))[0][0]])(getattr('whfg '
+'''nabgure clguba unpxre''', dir('0o52')[0o107])({ _:
(_-0o124) %0o32 +0o141 for _ in range(0o141, 0o173)}))

reformed japhs in python: ridiculous anagram

2011-04-03T10:30:00.000-04:00

Here's the second in my reformed JAPH series. It takes an anagram of 'just another python hacker' and converts it prior to printing. It sorts the anagram by the indices of another string, in order of their associated characters. This is sort of like a pre-digested Schwartzian transform.

x = 'upjohn tehran hectors katy'
y = '1D0HG6JFO9P5ICKAM87B24NL3E'

print(''.join(x[i] for i in sorted(range(len(x)), key=lambda p: y[p])))

Obfuscation consists mostly of using silly machinations to construct the string we use to sort the anagram.

print(''.join('''upjohn tehran hectors katy'''[_]for _ in sorted(range
(26),key=lambda p:(hex(29)[2:].upper()+str(3*3*3*3-3**4)+'HG'+str(sum(
range(4)))+'JFO'+str((1+2)**(1+1))+'P'+str(35/7)[:1]+'i.c.k.'.replace(
'.','').upper()+'AM'+str(3**2*sum(range(5))-3)+hex(0o5444)[2:].replace
(*'\x62|\x42'.split('|'))+'NL'+hex(0o076).split('x')[1].upper())[p])))

reformed japhs in python: alphabetic indexing

2011-04-01T21:40:00.008-04:00

Full disclosure: I used to be a Perl programmer. (Hey, Jeff.)

One day I became disillusioned at the progress of Perl 6 and decided to import this. This appears to be a fairly common story for Perl to Python converts. While I haven't looked back much, there are a number of things I really miss about perl (lower case intentional). Among other things, I miss having value types in a dynamic language, magical and ill-advised use of cryptocontext, and sometimes even pseudohashes because they were inexcusably weird. A language that supports so many ideas out of the box enables an extended learning curve that lasts for many years. "Perl itself is the game."

Most of all I think I miss writing Perl poetry and JAPHs. Stupidly, I didn't keep any of those I wrote, and I'm not competent enough with the language anymore to write interesting ones. At the time I was intentionally distancing myself from a model that was largely implicit and based on understanding of archaic system internals and moving to one that was (supposedly) explicit and simple. After switching to Python as my primary language, I used the following email signature in a nod to this change in worldview (intended for Python 2.x):

print 'just another python hacker'

Recently I've been experimenting with writing JAPHs in Python. I think of these as reformed JAPHs. They accomplish the same purpose as programming exercises but exist with a far more restricted context. In some ways they are more challenging. Creativity is difficult when functioning within a more explicit landscape. I now have a small series of reformed JAPHs dedicated to a close friend which increase monotonically in complexity. Here is the first one, written in plain understandable Python 3.2:

import string

letters = string.ascii_lowercase + ' '
indices = [
     9, 20, 18, 19, 26,  0, 13, 14, 19, 7,  4, 17, 26,
    15, 24, 19,  7, 14, 13, 26,  7,  0, 2, 10,  4, 17
]

print(''.join(letters[i] for i in indices))

This is fairly simple. Instead of explicitly embedding the string 'just another python hacker' in the program, we assemble it using the index of its letters in the string 'abcdefghijklmnopqrstuvwxyz '. Obfuscation is achieved through a series of minor measures:

Instead of calling the print function, we import sys and make a call to sys.stdout.write.
string.lowercase + ' ' is assembled by joining together the character versions of its respective ordinal values (97 to 123 and 32).
The integer indices are joined together using 'l' and split into a list.
Liberal use of ''' and the fact that multiple subsequent strings in Python are concatenated.

Here's the obfuscated version:

eval("__import__('''\x73''''''\x79''''''\x73''').sTdOuT".lower()
).write(''.join(map(lambda _:(list(map(chr,range(97,123)))+[chr(
32)])[int(_)],('''9l20l18l19''''''l26l0l13l14l19l7l4l17l26l15'''
'''l24l19l7l14l1''''''3l26l7l0l2l10l4l17''').split('l')))+'\n',)

More could certainly be done, but that's about as far as this one captivated my interest. Stay tuned for the next JAPH.

what i hate about writing decision support systems: a call to arms

2011-03-14T17:28:00.000-04:00

I realize it's been a while since I've updated python-zibopt. While it's true that I am wracked with guilt over this¹, I can't erase a feeling that there are larger issues that should be addressed in the optimization solver API department.

There are two standard methods for accessing a solver: through an algebraic modeling language like AMPL or using some kind of API. python-zibopt belongs to the latter category, along with a host of other technologies. Modeling languages tend to be very good at describing algebraic models succinctly, and at loading data into these models (provided it's in a supported format). Where they lack is in basic programming features and libraries. Try writing an interface between a modeling system and a web site, or performing anything beyond the simplest logic within a model. These are often difficult or impossible tasks. Never mind getting access to things like an object system or concurrency.

And these aren't things modeling language developers should have to worry about. They already exist in most general purpose programming languages, which is why most serious decision support systems use direct APIs for optimization. Unfortunately, this means giving up access to most algebraic modeling features. Consider the flow of operations for a typical decision support system:

Collect user input, probably represented in a database model.
Convert this input to an algebraic model and variables.
Hand this model off to a solver, often in the form of a big, ugly matrix.
Call the solver. Wait for an optimal solution, maximum bounds gap, or time out.
Assuming a feasible or optimal solution is even found, convert it back.
Display model solution output to the user.

This is a little bit insane. At best, a developer has to code two representations of the data and two different conversions from one model type to another. The pain of this process has a lot to do with the solver APIs available, and these are anything but standard. At its most painful (and most common), this requires construction of a giant matrix for feeding into the solver. Some APIs wrap constraints and "variables"² into object classes to try and attenuate the pain. Others, like python-zibopt, add syntactic sugar from algebraic modeling languages, such as creating constraints by summation of decision variables over index sets, to make the modeler feel more at home.

Still, that's about as far as we can seem to go. And the results still require different representations of the same information along with numerous conversions back and forth between formats. Surely, as operations research practitioners we recognize that this is fundamentally inefficient. As someone who straddles the worlds of CS and OR, here are a few things that I think should be easy for me to do when developing models:

Define single variables that represent data living in a database, web form, and solver, much the way web programmers use an object-relational mapper.
Serialize an algebraic model to a simple format such as JSON.
Subject decision variables to the same kinds of conditionals I can use with normal variables, like "if i <= 3".
All the things I haven't thought of yet but can do using a traditional programming language.

It's my belief that this is difficult because general purpose languages and modeling languages treat variables in fundamentally different ways. In most languages, the moment a variable is defined it has some sort of value, even if that value is null. In purely functional languages, variables are prohibited from even being null. (They also don't vary, but that's another story.) For an algebraic model, a variable represents a universe of possibilities. Typically a variable can take on any real value immediately upon instantiation. That universe may be limited by the addition of constraints. Decision variables then take on a series of values as the solver searches the resulting feasible set.

It seems that resolving this dichotomy could alleviate the difficulties of our current modeling situation. I'm not sure where to start. Perhaps the introduction of a sub-language into a general purpose language would do the trick? This would be similar to how regular expressions are treated in Perl 6. I'd like to hear your thoughts.

Seriously. This keeps me up at night.

The oddity of defining another variable type within a general purpose programming language is not lost on me.

talks that shaped me

2011-03-12T23:09:00.000-05:00

As this weekend's PyCon progresses, which I am unfortunately unable to attend, I find myself reminiscing about a similar conference I went to a decade ago. At the time the .com bust hadn't entirely taken hold of the Washington economy yet, so a few coworkers and I were sent to YAPC 19101¹. I was a newly minted developer and hadn't received much formal computer science education yet. Here are a few talks that opened my eyes to the world of possibilities.

Michael Schwern: Disciplined Programming, or How to be Lazy without Really Trying

This talk was about testing. Schwern attempted to convince us that writing unit tests, although they weren't called that yet, was fundamental to the three virtues of a programmer. He reasoned that he maintained a lot of well used modules on CPAN and would never be able to handle them all if he hadn't developed habits that enable laziness. Thus automated testing engenders laziness and is therefore to be encouragd. This is when I realized that programming is the intersection of two totally different fields: a science and a discipline.

Mark-Jason Dominus: Stolen Secrets of the Wizards of the Ivory Tower

MJD, one of the few programmers known by his initials, was deep into his studies of functional programming at the time and gave an in-depth preview of what would become Higher Order Perl. In addition to well known advice from ESR, another TLA, this is why I learned functional programming. It's also what piqued my interest in algorithmic complexity, which later led to a love of integer programming.

Damian Conway: Life, the Universe, and Everything

This talk is mostly a blur, but that's because it was at the time too. Conway possessed such a deep understanding of so many topics, and an ability to leap from one to another, that the experience was dizzying. He implemented the Game of Life using Quantum::Superpositions, both wrote Perl and spoke in Klingon, and, if I remember correctly, received a standing ovation at the end. Conway set my standard for live demos.

There were others that made an impression, but even these three are hazy after a decade of bit rot.

When I got home I was excited and in awe. I hacked at tail recursive functions in assembly language, studied LISP, and painstakingly implemented the exercises from SICP in half a dozen languages. I learned Erlang from the only available book at the time because I thought concurrent programming would one day be important². I picked up Python because "Python is not the enemy, Java is." YAPC energized me to pursue the routine of computer science autodidact. It took almost another ten years.

Montreal will always hold a special place in my heart. Before I was born, my parents lived there while my father was a graduate student at McGill studying biophysical chemistry. They were dirt poor and suffered, but whenever they speak of it they light up and are suddenly animated. For me, it's where I learned how big the universe of knowledge is, and how little I had been exposed to. It's where I was given the push that shaped my studies, career, and life. Nowadays conference talks are often recorded and posted to YouTube or blip.tv, but back then it was considered cutting edge for someone to post MP3s of lighting talks. To those of you who speak at community technical conferences, thank you, and know that your efforts do make an impact.

Prior to Y2K, it was fairly common to find Perl code that assembled the current year using "'19'.(localtime())[5]". Of course, index 5 of the return value of localtime is the number of years since 1900, so this produced 19100 instead of 2000 and 19101 instead of 2001. This was really funny at the time.

Though it's uncredited, I actually did submit an entry for OTP to FOLDOC back in August of 2001. My claim to fame.

simulating gdp growth

2011-02-23T17:11:00.004-05:00

[2011-02-24 This post has been corrected. Near the end it stated that we could be 95% confident China's GDP would overtake the USA between 2057 and 2059. It should have stated that we could be 95% confident the mean of that year would be between 2057 and 2059. Thanks to my overbearing subconscious for pointing this out.]

I hope you saw China’s way to the top on the Post's website recently. It's a very clear presentation of their statement and is certainly worth a look.

So say you're an economist and you actually do need to produce a realistic estimate of when China's GDP surpasses that of the USA. Can you use such an approach? Not really. There are several simplifying assumptions the Post made that are perfectly reasonable. However, if the goal is an analytical output from a highly random system such as GDP growth, one should not assume the inputs are fixed. (I'm not saying I have any gripe with their interactive. This post has a different purpose.)

Why is this? The short answer is that randomness in any system can change its output drastically from one run to the next. Even if the mean from a deterministic analysis is correct, it tells us nothing about the variance of our output. We really need a confidence interval of years when China is likely to overtake the USA.

We'll move in the great tradition of all simulation studies. First we pepare our input. A CSV of GDP in current US dollars for both countries from 1960 to 2009 is available here. This was simply copied and transposed from the World Bank data files (thanks, Kelso). We read this into a data frame and calculate their growth rates year over year. Note that the first value for growth has to be NA.

> gdp <- read.csv('gdp.csv')
> gdp$USA.growth <- rep(NA, length(gdp$USA))
> gdp$China.growth <- rep(NA, length(gdp$China))
> for (i in 2:length(gdp$USA)) {
+     gdp$USA.growth[i] <- 100 * (gdp$USA[i] - gdp$USA[i-1]) / gdp$USA[i-1]
+     gdp$China.growth[i] <- 100 * (gdp$China[i] - gdp$China[i-1]) / gdp$China[i-1]
+ }

We now analyze our inputs and assign probability distributions to the annual growth rates. In a full study this would involve comparing a number of different distributions and choosing the one that fits the input data best, but that's well beyond the scope of this post. Instead, we'll use the poor man's way out: plot histograms and visually verify what we hope to be true, that the distributions are normal.

And they pretty much are. That's good enough for our purposes. Now all we need are the distribution parameters, which are mean and standard deviation for normal distributions.

> mean(gdp$USA.growth[!is.na(gdp$USA.growth)])
[1] 7.00594
> sd(gdp$USA.growth[!is.na(gdp$USA.growth)])
[1] 2.889808
> mean(gdp$China.growth[!is.na(gdp$China.growth)])
[1] 9.90896
> sd(gdp$China.growth[!is.na(gdp$China.growth)])
[1] 10.5712

Now our input analysis is done. These are the inputs:

USA Growth ~ N(7.00594, 2.889808²)
China Growth ~ N(9.90896, 10.5712²)

This should make the advantage of such an approach much more obvious. Compare the standard deviations for the two countries. China is a lot more likely to have negative GDP growth in any given year. They're also more likely to have astronomical growth.

We now build and run our simulation study. The more times we run the simulation the tighter we can make our confidence interval (to a point), so we'll pick a pretty big number somewhat arbitrarily. If we want to, we can be fairly scientific about determining how many iterations are necessary after we've done some runs, but we have to start somewhere.

> repetitions <- 10000

This is the code for our simulation. For each iteration, it starts both countries at their 2009 GDPs. It then iterates, changing GDP randomly until China's GDP is at least the same value as the USA's. When that happens, it records the current year.

> results <- rep(NA, repetitions)
> for (i in 1:repetitions) {
+     usa <- gdp$USA[length(gdp$USA)]
+     china <- gdp$China[length(gdp$China)]
+     year <- gdp$Year[length(gdp$Year)]
+  
+     while (TRUE) {
+         year <- year + 1
+   
+         usa.growth <- rnorm(1, 7.00594, 2.889808)
+         china.growth <- rnorm(1, 9.90896, 10.5712)
+   
+         usa <- usa * (1 + (usa.growth / 100))
+         china <- china * (1 + (china.growth / 100))
+   
+         if (china >= usa) {
+             results[i] <- year
+             break
+         }
+     }
+ }

From the results vector we see that, given the data and assumptions for this model, China should surpass the USA in 2058. We also see that we can be 95% confident that the mean year this will happen is between 2057 and 2059. This is not quite the same as saying we are confident this will actually happen between those years. The result of our simulation is a probability distribution and we are discovering information about it.

> mean(results)[1] 2058.494
> mean(results) + (sd(results) / sqrt(length(results)) * qnorm(0.025))
[1] 2057.873
> mean(results) + (sd(results) / sqrt(length(results)) * qnorm(0.975))
[1] 2059.114

So what's wrong with this model? Well, we had to make a number of assumptions:

We assume we actually used the right data set. This was more of a how-to than a proper analysis, so that wasn't too much of a concern.
We assume future growth for the next 40-50 years resembles past growth from 1960-2009. This is a bit ridiculous, of course, but that's the problem with forecasting.
We assume growth is normally distributed and that we don't encounter heavy-tailed behaviors in our distributions.
We assume each year's growth is independent of the year before it. See the last exercise.

Here are some good simulation exercises if you're looking to do more:

Note how the outputs are quite a bit different from the Post graphic. I expect that's largely due to the inclusion of data back to 1960. Try running the simulation for yourself using just the past 10, 20, and 30 years and see how that changes the result.
Write a simulation to determine the probability China's GDP surpasses the USA's in the next 25 years. Now plot the mean GDP and 95% confidence intervals for each country per year.
Assume that there are actually two distributions for growth for each country: one when the previous year had positive growth and another when it was negative. How does that change the output?

data fitting part 2a: very, very simple linear regression in r

2011-02-15T16:13:00.000-05:00

I thought it might be useful to follow up the last post with another one showing the same examples in R. If you want to try the examples, first download a CSV of the input data here.

R provides a function called lm, which is similar in spirit to numpy's linalg.lstsq. As you'll see, lm's interface is a bit more tuned to the concepts of modeling.

We begin by reading in the example CSV into a data frame:

> responses <- read.csv('example_data.csv')

> responses
  respondent vanilla.love strawberry.love chocolate.love
1     Serdar            9               4              9
2        Dan            8               6              4
3  Nathaniel            9               4              8
4     Lauren            3               7              9
5        Jen            6               8              5
6     Jackie            4               5              3

  dog.love cat.love
1        9        9
2       10        4
3        2        6
4        4        6
5        2        5
6       10        3

A data frame is sort of like a matrix, but with named columns. That is, we can refer to entire columns using the dollar sign. We are now ready to run least squares. We'll create the model for predicting "dog love." To create the "cat love" model, simply use that column name instead:

> fit1 <- lm(responses$dog.love ~ responses$vanilla.love 
    + responses$strawberry.love + responses$chocolate.love)

The syntax for lm is a little offputting at first. This call tells it to create a model for "dog love" with respect to (the ~) a function of the form offset + x1*vanilla love + x2*strawberry love + x3*chocolate love. Note that the offset is conveniently implied when using lm, so this is the same as the second model we created in Python. Now that we've computed the coefficients for our "dog love" model, we can ask R about it:

> summary(fit1)

Call:
lm(formula = responses$dog.love ~ responses$vanilla.love 
  + responses$strawberry.love + responses$chocolate.love)

Residuals:
      1       2       3       4       5       6 
 3.1827  2.9436 -4.5820  0.8069 -1.9856 -0.3657 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)
(Intercept)                20.9298    15.0654   1.389    0.299
responses$vanilla.love     -0.2783     0.9934  -0.280    0.806
responses$strawberry.love  -1.4314     1.5905  -0.900    0.463
responses$chocolate.love   -0.7647     0.8214  -0.931    0.450

Residual standard error: 4.718 on 2 degrees of freedom
Multiple R-squared: 0.4206,     Adjusted R-squared: -0.4485 
F-statistic: 0.484 on 3 and 2 DF,  p-value: 0.7272

This gives us quite a bit of information, including the coefficients for our "dog love" model and various error metrics. You can find the offset and coefficients under the Estimate column above. We quickly verify this using R's vectorized arithmetic:

> 20.9298 - 0.2783 * responses$vanilla.love 
    - 1.4314 * responses$strawberry.love 
    - 0.7647 * responses$chocolate.love

[1]  5.8172  7.0562  6.5819  3.1928  3.9853 10.3655

You'll notice the model is essentially the same as the one we got from numpy. Our next step is to add in the squared inputs. We do this by adding extra terms to the modeling formula. The I() function allows us to easily add additional operators to columns. That's how we accomplish the squaring. We could alternatively add squared input values to the data frame, but using I() is more convenient and natural.

> fit2 <- lm(responses$dog.love ~ responses$vanilla.love 
    + I(responses$vanilla.love^2) + responses$strawberry.love 
    + I(responses$strawberry.love^2) + responses$chocolate.love 
    + I(responses$chocolate.love^2))

> summary(fit2)

Call:
lm(formula = responses$dog.love ~ responses$vanilla.love 
  + I(responses$vanilla.love^2) + responses$strawberry.love 
  + I(responses$strawberry.love^2) + responses$chocolate.love 
  + I(responses$chocolate.love^2))

Residuals:
ALL 6 residuals are 0: no residual degrees of freedom!

Coefficients: (1 not defined because of singularities)
                               Estimate Std. Error t value Pr(>|t|)
(Intercept)                    -357.444         NA      NA       NA
responses$vanilla.love           72.444         NA      NA       NA
I(responses$vanilla.love^2)      -6.111         NA      NA       NA
responses$strawberry.love        59.500         NA      NA       NA
I(responses$strawberry.love^2)   -5.722         NA      NA       NA
responses$chocolate.love          7.000         NA      NA       NA
I(responses$chocolate.love^2)        NA         NA      NA       NA

Residual standard error: NaN on 0 degrees of freedom
Multiple R-squared:     1,      Adjusted R-squared:   NaN 
F-statistic:   NaN on 5 and 0 DF,  p-value: NA

We can see that we get the same "dog love" model as produced by the third Python version of the last post. Again, we quickly verify that the output is the same (minus some rounding errors):

> -357.444 + 72.444 * responses$vanilla.love 
    - 6.111 * responses$vanilla.love^2 
    + 59.5 * responses$strawberry.love 
    - 5.722 * responses$strawberry.love^2 
    + 7 * responses$chocolate.love

[1]  9.009 10.012  2.009  4.011  2.016 10.006

data fitting part 2: very, very simple linear regression in python

2011-02-15T00:39:00.012-05:00

This post is based on a memo I sent to some former colleagues at the Post. I've edited it for use here since it fits well as the second in a series on simple data fitting techniques. If you're among the many enlightened individuals already using regression analysis, then this post is probably not for you. If you aren't, then hopefully this provides everything you need to develop rudimentary predictive models that yield surprising levels of accuracy.

For purposes of a simple working example, we have collected six records of input data over three dimensions with the goal of predicting two outputs. The input data are:


    x1: How much a respondent likes vanilla [0-10]
    x2: How much a respondent likes strawberry [0-10]
    x3: How much a respondent likes chocolate [0-10]

Output data consist of:


    b1: How much a respondent likes dogs [0-10]
    b2: How much a respondent likes cats [0-10]

Below are entirely anonymous data collected from people who bear no resemblance to certain Washington Post staffers.

respondent	vanilla love	strawberry love	chocolate love	dog love	cat love
Serdar	9	4	9	9	9
Dan	8	6	4	10	4
Nathaniel	9	4	8	2	6
Lauren	3	7	9	4	6
Jen	6	8	5	2	5
Jackie	4	5	3	10	3

Our input is in three dimensions. Each output requires its own model, so we'll have one for dogs and one for cats. We're looking for functions, dog(x) and cat(x), that can predict b1 and b2 based on given values of x1, x2, and x3.

For both models we want to find parameters that minimize their squared residuals (read: errors). There's a number of names for this. Optimization folks like to think of it as unconstrained quadratic optimization, but it's more common to call it least squares or linear regression. It's not necessary to entirely understand why for our purposes, but the function that minimizes these errors is:

$\beta = ({A^t}A)^{-1}{A^t}b$

This is implemented for you in the numpy.linalg Python package, which we'll use for examples. Much more information than you probably want can be found here.

Below is a first stab at a Python version. It runs least squares against our input and output data exactly as they are. You can see the matrix A and outputs b1 and b2 (dog and cat love, respectively) are represented just as they are in the table.

# Version 1: No offset, no squared inputs

import numpy

A = numpy.vstack([
    [9, 4, 9],
    [8, 6, 4],
    [9, 4, 8],
    [3, 7, 9],
    [6, 8, 5],
    [4, 5, 3]
])

b1 = numpy.array([9, 10, 2, 4, 2, 10])
b2 = numpy.array([9, 4, 6, 6, 5, 3])

print 'dog <3:', numpy.linalg.lstsq(A, b1)[0]
print 'cat <3:', numpy.linalg.lstsq(A, b2)[0]

# Output:
# dog <3: [0.72548294      0.53045642     -0.29952361]
# cat <3: [2.36110929e-01  2.61934385e-05  6.26892476e-01]

The resulting model is:

dog(x) = 0.72548294 * x1 + 0.53045642 * x2 - 0.29952361 * x3
cat(x) = 2.36110929e-01 * x1 + 2.61934385e-05 * x2 + 6.26892476e-01 * x3

The coefficients before our variables correspond to beta in the formula above. Errors between observed and predicted data, shown below, are calculated and summed. For these six records, dog(x) has a total error of 20.76 and cat(x) has 3.74. Not great.

Person	predicted b1	b1 error	predicted b2	b2 error
Serdar	5.96	3.04	7.77	1.23
Dan	7.79	2.21	4.40	0.40
Nathaniel	6.25	4.25	7.14	1.14
Lauren	3.19	0.81	6.35	0.35
Jen	7.10	5.10	4.55	0.45
Jackie	4.66	5.34	2.83	0.17
Total error:		20.76		3.74

One problem with this model is that dog(x) and cat(x) are forced to pass through the origin. (Why is that?) We can improve it somewhat if we add an offset. This amounts to prepending 1 to every row in A and adding a constant to the resulting functions. You can see the very slight difference between the code for this model and that of the previous:

# Version 2: Offset, no squared inputs

import numpy

A = numpy.vstack([
    [1, 9, 4, 9],
    [1, 8, 6, 4],
    [1, 9, 4, 8],
    [1, 3, 7, 9],
    [1, 6, 8, 5],
    [1, 4, 5, 3]
])

b1 = numpy.array([9, 10, 2, 4, 2, 10])
b2 = numpy.array([9, 4, 6, 6, 5, 3])

print 'dog <3:', numpy.linalg.lstsq(A, b1)[0]
print 'cat <3:', numpy.linalg.lstsq(A, b2)[0]

# Output:
# dog <3: [20.92975427  -0.27831197  -1.43135684  -0.76469017]
# cat <3: [-0.31744124   0.25133547   0.02978098   0.63394765]

This yields the seconds version of our models:

dog(x) = 20.92975427 - 0.27831197 * x1 - 1.43135684 * x2 - 0.76469017 * x3
cat(x) = -0.31744124 + 0.25133547 * x1 + 0.02978098 * x2 + 0.63394765 * x3

These models provide errors of 13.87 and 3.79. A little better on the dog side, but still not quite usable.

Person	predicted b1	b1 error	predicted b2	b2 error
Serdar	5.82	3.18	7.77	1.23
Dan	7.06	2.94	4.41	0.41
Nathaniel	6.58	4.58	7.14	1.14
Lauren	3.19	0.81	6.35	0.35
Jen	3.99	1.99	4.60	0.40
Jackie	10.37	0.37	2.74	0.26
Total error:		13.87		3.79

The problem is that dog(x) and cat(x) are linear functions. Most observed data don't conform to straight lines. Take a moment and draw the line $f(x) = x$ and the curve $f(x) = x^2$. The former makes a poor approximation of the latter.

Most of the time, people just use squares of the input data to add curvature to their models. We do this in our next version of the code by just adding squares of the input row values to our A matrix. Everything else is the same. (In reality, you can add any function of the input data you feel best models the data, if you understand it well enough.)

# Version 3: Offset with squared inputs

import numpy

A = numpy.vstack([
    [1, 9, 9**2, 4, 4**2, 9, 9**2],
    [1, 8, 8**2, 6, 6**2, 4, 4**2],
    [1, 9, 9**2, 4, 4**2, 8, 8**2],
    [1, 3, 3**2, 7, 7**2, 9, 9**2],
    [1, 6, 6**2, 8, 8**2, 5, 5**2],
    [1, 4, 4**2, 5, 5**2, 3, 3**2]
])

b1 = numpy.array([9, 10, 2, 4, 2, 10])
b2 = numpy.array([9, 4, 6, 6, 5, 3])

print 'dog <3:', numpy.linalg.lstsq(A, b1)[0]
print 'cat <3:', numpy.linalg.lstsq(A, b2)[0]

# dog <3: [1.29368307  7.03633306  -0.44795498  9.98093332  
#  -0.75689575  -19.00757486  1.52985734]
# cat <3: [0.47945896  5.30866067  -0.39644128 -1.28704188  
#   0.12634295   -4.32392606  0.43081918]

This gives us our final version of the model:

dog(x) = 1.29368307 + 7.03633306 * x1 - 0.44795498 * x1**2 + 9.98093332 * x2 - 0.75689575 * x2**2 - 19.00757486 * x3 + 1.52985734 * x3**2
cat(x) = 0.47945896 + 5.30866067 * x1 - 0.39644128 * x1**2 - 1.28704188 * x2 + 0.12634295 * x2**2 - 4.32392606 * x3 + 0.43081918 * x3**2

Adding curvature to our model eliminates all perceived error, at least within 1e-16. This may seem unbelievable, but when you consider that we only have six input records, it isn't really.

Person	predicted b1	predicted b2
Serdar	9	9
Dan	10	4
Nathaniel	2	6
Lauren	4	6
Jen	2	5
Jackie	10	3
Total error:

It should be fairly obvious how one can take this and extrapolate to much larger models. I hope this is useful and that least squares becomes an important part of your lives.

off-the-cuff online voter fraud detection

2010-11-30T23:36:00.000-05:00

Consider this scenario: You run a contest that accepts votes from the general Internet population. In order to encourage user engagement, you record any and all votes into a database over several days, storing nothing more than the competitor voted for, when each vote is cast, and a cookie set on the voter's computer along with their apparent IP addresses. If a voter already has a recorded cookie set they are denied subsequent votes. This way you can avoid requiring site registration, a huge turnoff for your users. Simple enough.

Unfortunately, some of the competitors are wily and attached to the idea of winning. They go so far as programming or hiring bots to cast thousands of votes for them. Your manager wants to know which votes are real and which ones are fake Right Now. Given very limited time, and ignoring actions that you could have taken to avoid the problem, how can you tell apart sets of good votes from those that shouldn't be counted?

One quick-and-dirty option involves comparing histograms of interarrival times for sets of votes. Say you're concerned that all the votes during a particular period of time or from a given IP address might be fraudulent. Put all the vote times you're concerned about into a list, sort them, and compute their differences:

# times is a list of datetime instances from vote records
times.sort(reversed=True)
interarrivals = [y-x for x, y in zip(times, times[1:]]

Now use matplotlib to display a histogram of these. Votes that occur naturally are likely to resemble an exponential distribution in their interarrival times. For instance, here are interarrival times for all votes received in a contest:


Interarrival times for all submissions

This subset of votes is clearly fraudulent, due to the near determinism of their interarrival times. This is most likely caused by the voting bot not taking random sleep intervals during voting. It casts a vote, receives a response, clears its cookies, and repeats:

Interarrival times for clearly fraudulent votes

These votes, on the other hand, are most likely legitimate. They exhibit a nice Erlang shape and appear to have natural interarrival times that one would expect:

Proper-looking interarrival times

Of course this method is woefully inadequate for rigorous detection of voting fraud. Ideally one would find a method to compute the probability that a set of votes is generated by a bot. This is enough to inform quick, ad hoc decisions though.

data fitting part 1: linear data fitting

2010-11-23T18:00:00.002-05:00

Data fitting is one of those tasks that everyone should have at least some exposure to. Certainly developers and analysts will benefit from a working knowledge of its fundamentals and their implementations. However, in my own reading I've found it difficult to locate good examples that are simple enough to pick up quickly and come with accompanying source code.

This article commences an ongoing series introducing basic data fitting techniques. With any luck they won't be overly complex, while still being useful enough to get the point across with a real example and real data. We'll start with a binary classification problem: presented with a series of records, each containing a set number of input values describing it, determine whether or not each record exhibits some property.

We'll use the cancer1.dt data from the proben1 set of test cases, which you can download here. Each record starts with 9 data points containing physical characteristics of a tumor. The second to last data point contains 1 if a tumor is benign and 0 if it is malignant. We seek to find a linear function we can run on an arbitrary record that will return a value greater than zero if that record's tumor is predicted to be benign and less than zero if it is predicted to be malignant. We will train our linear model on the first 350 records, and test it for accuracy on the remaining rows.

This is similar to the data fitting problem found in Chvatal. Our inputs consist of a matrix of observed data, $A$, and a vector of classifications, $b$. In order to classify a record, we require another vector $x$ such that the dot product of $x$ and that record will be either greater or less than zero depending on its predicted classification.

A couple points to note before we start:

Most observed data are noisy. This means it may be impossible to locate a hyperplane that cleanly separates given records of one type from another. In this case, we must resort to finding a function that minimizes our predictive error. For the purposes of this example, we'll minimize the sum of the absolute values of the observed and predicted values. That is, we seek $x$ such that we find $min \sum_i{|a_i^T x-b_i|}$.
The slope-intercept form of a line, $f(x)=m^T x+b$, contains an offset. It should be obvious that this is necessary in our model so that our function isn't required to pass through the origin. Thus, we'll be adding an extra variable with the coefficient of 1 to represent our offset value.

In order to model this, we use two linear constraints for each absolute value. We minimize the sum of these. Our Linear Programming model thus looks like:

$\min$	$z = x_0 + \sum_i{v_i}$
s.t.	$v_i$	$\geq$	$x_0 + a_i'x - 1$	$\forall$ benign tumors
	$v_i$	$\geq$	$1 - x_0 - a_i'x$	$\forall$ benign tumors
	$v_i$	$\geq$	$x_0 + a_i'x - (-1)$	$\forall$ malignant tumors
	$v_i$	$\geq$	$-1 - x_0 - a_i'x$	$\forall$ malignant tumors

In order to do this in Python, we use SCIP and SoPlex from the zibopt library. We start by setting constants for benign and malignant outputs and providing a function to read in the training and testing data sets.

# Preferred output values for tumor categories
BENIGN = 1
MALIGNANT = -1

def read_proben1_cancer_data(filename, train_size):
    '''Loads a proben1 cancer file into train & test sets'''
    # Number of input data points per record
    DATA_POINTS = 9

    train_data = []
    test_data = []

    with open(filename) as infile:
        # Read in the first train_size lines to a training
        # data list, and the others to testing data.  This
        # allows us to test how general our model is on
        # something other than the input data.
        for line in infile.readlines()[7:]: # skip header
            line = line.split()

            # Records = offset (x0) + remaining data points
            input = [float(x) for x in line[:DATA_POINTS]]
            output = BENIGN if line[-2] == '1' else MALIGNANT
            record = {'input': input, 'output': output}

            # Determine what data set to put this in
            if len(train_data) >= train_size:
                test_data.append(record)
            else:
                train_data.append(record)

    return train_data, test_data

The next function implements the LP model described above using SoPlex and SCIP. It minimizes the sum of residuals for each training record. This amounts to summing the absolute value of the difference between predicted and observed output data. The following function takes in input and observed output data and returns a list of coefficients. Our resulting model consists of taking the dot product of an input record and these coefficients. If the result is greater than or equal to zero, that record is predicted to be a benign tumor, otherwise it is predicted to be malignant.

def train_linear_model(train_data):
    '''
    Accepts a set of input training data with known output
    values.  Returns a list of coefficients to apply to
    arbitrary records for purposes of binary categorization.
    '''
    from zibopt import scip
    import sys

    # Make sure we have at least one training records
    assert len(train_data) > 0
    num_variables = len(train_data[0]['input'])

    # Variables are coefficients in front of the data points.
    # It is important that these be unrestricted in sign so
    # they can take negative values.
    solver = scip.solver()
    coefficients = [
        solver.variable(lower=-sys.maxint)
        for _ in xrange(num_variables)
    ]

    # Residual for each data row
    residuals = [solver.variable() for _ in train_data]
    for r, d in zip(residuals, train_data):
        # r will be the absolute value of the difference
        # between observed and predicted values.  We can
        # model absolute values such as r >= |foo| as:
        #
        #   r >=  foo
        #   r >= -foo
        solver += sum(
            x * y for x, y in zip(coefficients, d['input'])
        ) + r >= d['output']

        solver += sum(
            x * y for x, y in zip(coefficients, d['input'])
        ) - r <= d['output']

    # Find and return coefficients that min sum of residuals
    solution = solver.minimize(objective=sum(residuals))
    return [solution[c] for c in coefficients]

We also provide a convenience function for counting the number of correct predictions by our resulting model against either the test or training data sets.

def count_correct(data_set, coefficients):
    '''Returns the number of correct predictions'''
    correct = 0
    for d in data_set:
        result = sum(
            x*y for x, y in zip(coefficients, d['input'])
        )

        # Do we predict the same as the output?
        if (result >= 0) == (d['output'] >= 0):
            correct += 1

    return correct

Finally we write a main method to read in the data, build our linear model, and test its efficacy.

if __name__ == '__main__':
    from pprint import pprint

    # Specs for this input file
    INPUT_FILE_NAME = 'cancer1.dt'
    TRAIN_SIZE = 350

    train_data, test_data = read_proben1_cancer_data(
        INPUT_FILE_NAME, TRAIN_SIZE
    )

    # Add the offset variable to each of our data records
    for data_set in [train_data, test_data]:
        for row in data_set:
            row['input'] = [1] + row['input']

    coefficients = train_linear_model(train_data)
    print 'coefficients:'
    pprint(coefficients)

    # Print % of correct preditions for each data set
    correct = count_correct(train_data, coefficients)
    print '%s / %s = %.02f%% correct on training set' % (
        correct, len(train_data),
        100 * float(correct) / len(train_data)
    )

    correct = count_correct(test_data, coefficients)
    print '%s / %s = %.02f%% correct on testing set' % (
        correct, len(test_data),
        100 * float(correct) / len(test_data)
    )

The result of running this model against the cancer1.dt data set is:

coefficients:
[1.4072882449702786,
 -0.14014055927954652,
 -0.6239513714263405,
 -0.26727681774258882,
 0.067107753841131157,
 -0.28300216102808429,
 -1.0355594670918404,
 -0.22774451038152174,
 -0.69871243677663608,
 -0.072575089848659444]
328 / 350 = 93.71% correct on training set
336 / 349 = 96.28% correct on testing set

The accuracy is pretty good here against the both the training and testing sets, so this particular model generalizes well. This is about the simplest model we can implement for data fitting, and we'll get to more complicated ones later, but it's nice to see we can do so well so quickly. The coefficients correspond to using a function of this form, rounding off to three decimal places:

$$f(x) = 1.407 - 0.140 x_1 - 0.624 x_2 - 0.267 x_3 + 0.067 x_4 - 0.283 x_5 - 1.037 x_6 - 0.228 x_7 - 0.699 x_8 - 0.073 x_9$$

moving on

2010-11-20T01:47:00.001-05:00

It's now official and public, so I may as well state it here: I'm leaving the Editorial Embedded Development team at the Washington Post for a position as a simulation & modeling engineer at MITRE. This represents a fundamental shift of career for me from software development to operations research, something I've hoped for since I took my first course in deterministic models at GMU. That said, leaving the Post remains a difficult step. In my tenure there, I've worked with many truly remarkable people on projects of obvious importance. I've weathered a difficult newsroom merger and been trusted with leadership of a blossoming development team. I've worked 32 hours straight to meet editorial deadlines and make news, won journalism awards by proxy, mentored aspiring developers, and been continually humbled by the intelligence, hard work, and dedication of my coworkers to accuracy, relevance, and objectivity. It is, therefore, not with a light heart that I take this step.

However, I know that it is the right path. MITRE is a great place to be for someone like me, and with any luck I'll have both time and encouragement to continue my contributions to open source optimization. As a first step, I'm picking this blog back up and finding time to post more interesting optimization tidbits.

recent work & graduation

2010-05-29T00:30:00.000-04:00

I spent the last term finishing my MSOR, so I haven't had much time to post since PyCon. It's been a busy term, and while I'm glad my efforts of the last few years have reached their conclusion, I'm already a little nostalgic for the classroom. Our department at George Mason has some really top notch faculty and the program is designed well. The fact that it requires MS students to balance between deterministic and stochastic methods is part of this. On my own I probably wouldn't have taken much on the stochastic side, but this way I ended up challenged and delighted by courses in simulation, stochastic processes, and queuing theory. I've met students from other institutions that don't require such a degree of breadth, and it makes me doubly happy to have studied where I did.

Here are a couple highlights of the past few months:

In the capstone project course we assembled onto teams and spent the term working for a project sponsor. Our team of five smart guys chose to work with Dr. Kuo-Chu Chang, whose interest was in evaluating financial engineering strategies against S&P 500 options. It was like opening a fire hose of information and required a huge effort through the entire term, but our result was a pretty solid system and analysis. I've posted a copy of the project web site here and the final report here. I was, of course, responsible for the software powering our deterministic simulation and GUI.

For queuing theory I had the opportunity to learn from a real expert in the field, which is always exciting. I've uploaded a short paper analyzing a queuing simulating of the Django installation we use at the Post here. It may be of interest to web programmers and administrators.

So what's next?

Well, first will come a new version of python-zibopt. After that I may spend a little time polishing PyGEP and then it's off to building our combinatorial auction software. Check back for updates!

my pycon picks for 2010

2010-02-09T15:34:00.001-05:00

Friday, Feb 19:

Saturday, Feb 20:

Sunday, Feb 21:

Dealing with unsightly data in the real world.
Optimal Resource Allocation using Python (of course I have to go to this one...)

pycon 2010 slides

2010-02-09T14:53:00.000-05:00

In case you'd like to see the working draft of my upcoming presentation on python-zibopt, you can find it here. I'm looking forward to seeing you there!

python-zibopt 0.4 released

2010-01-03T22:03:00.002-05:00

This is the last release I'll be making with new features before PyCon. That said, it's a pretty important one. First off, the clunky syntax for adding constraints has been replaced with a much more Pythonic one that allows algebraic notation using operator overloading:

from zibopt import scip
solver = scip.solver()

x1 = solver.variable(scip.INTEGER)
x2 = solver.variable(scip.INTEGER)
x3 = solver.variable(scip.INTEGER)

solver += x1 <= 2
solver += x1 + x2 + 3*x3 <= 3

solution = solver.maximize(objective=x1 + x2 + 2*x3)

It even allows you to do fancy things with sum(...) and to set both upper and lower bounds at the same time:

solver += 50 <= sum(cost[i] * x[i] for i in things) <= 100

The second new feature allows one to set branching priority on variables. This is pretty straightforward:

x1.priority = 1000

I'm still consulting with the ZIB folks about simplifying the installation, but I'm feeling pretty good about the state of the API. Time to start practicing my presentation!You can download the new release here.

my pycon 2010 talk: optimal resource allocation using python

2009-11-10T23:27:00.000-05:00

I'm very excited to report that my PyCon 2010 talk has been accepted and I will be publicly presenting python-zibopt in Atlanta this coming February. There's still some work to be done on the project, including simplifying the addition of constraints and figuring out an easy way to distribute it, but for now I give you the outline to:

Optimal Resource Allocation using Python

A brief introduction to modeling and solving resource allocation and scheduling problems using Python and SCIP.

At times a programmer is faced with difficult, possibly NP-Hard, optimization problems such as scheduling or assignment. Sophisticated techniques exist for modeling and solving these sorts of problems which are well implemented in optimization solvers. This talk introduces some of these techniques using the ZIB Optimization Suite and its new Python interface.

Outline

What is Combinatorial Optimization?
A few methods for finding solutions
OK, but how do I really do it?
Introducing python-zibopt
Conclusion

What is Combinatorial Optimization?

A brief history of Mathematical Optimization
Example of Linear Optimization: Production planning
- Problem statement & data
- Mathematical model
The Combinatorial Explosion: what happens when I can't produce 1/2 a widget?
- Example: Uncapacitated Facility Location (UFL)
- Example: The Traveling Salesman Problem
- What makes these problems difficult?

A few methods for finding solutions

It's all about bounding
- Primal bounds: I have a feasible solution and here is its objective value
- Dual bounds: I've relaxed my problem and can do no better than this
- Use Primal/Dual gap for search tree pruning
Branch-and-Bound
Prosolving (or preprocessing) on binary integer programs
Cutting planes: cut off fractional solutions from LP dual

OK, but how do I really do it?

You can always code it up yourself, but...
All of these techniques (and more) are parts of commercial and open source solvers
Two common methods to interface with a solver:
- Modeling language example: ZIMPL
  - The quickest way to get started
  - Closely resembles mathematical description of problem
  - Separates data from model
  - Not a full programming language... It is difficult to:
    - Automatically change the model depending on output (TSP subtour elimination)
    - Combine data from disparate sources (CSV, RDBMS, arbitrarily structured text...)
    - Base parts of the model on logic or user input
- Solver APIs
  - Give all the existing power of a solver
  - Allow you to build your own problem-specific optimization components
  - Are usually in C...

Introducing python-zibopt

A brief intro to the ZIB Optimization Suite
- Modeling Language: ZIMPL
- Linear Programming Solver: SoPlex
- Integer Programming Solver: SCIP
Features of python-zibopt
- Creating model variables
- Describing problem constraints
- Maximization / minimization
- Maximum gaps and time limits on solving
- Example: Maximizing network throughput
  - You have some kind of network (computer network, transportation, etc). Arcs have given capacities and you want to (a) calculate maximum throughput and (b) know where best to build extensions to the network.
- Example: Scheduling of on-call times for IT admins
  - You have servers and services running in an IT shop. Services go down in the middle of the night and have to be fixed and brought back up. You have limited IT staff, each with different capabilities (apache, oracle, python, etc). Given approved vacations, you want to create a feasible on-call schedule that equitably assigns on-call duties. You would also like guidance for future hiring and training decisions.

Conclusion

Just because a problem is intractable doesn't mean you can't solve it
ZIB has an excellent suite of open source tools for just this purpose
Optimization is for Python programmers too!

python-zibopt 0.2 released

2009-10-17T01:14:00.001-04:00

The new release has some exciting features, most importantly:

You can now use gap, absgap, and time keywords when calling maximize/minimize to stop the solver before reaching optimality. These function the same way they function for SCIP, pausing the search when any one of the conditions is met. You can use any combination of them, and on the next call to maximize or minimize SCIP will pick up where it left off.

solver.maximize(absgap=0.5, gap=100, time=500)

Settings for branching rules, separators, and conflict handlers can be read and modified from within Python:

solver.branching['inference'].priority = 10000
solver.branching['mostinf'].maxdepth += 2
print(solver.branching['pscost'].maxbounddist)

solver.separators['clique'].priority = 10000
solver.conflict['logicor'].priority = 10000

You can download the new release here.

easy monte carlo simulation in python

2009-10-08T12:11:00.008-04:00

One of the most useful tools one learns in an Operations Research curriculum is
Monte Carlo Simulation. Its utility lies in its simplicity: one can learn vital information about nearly any process, be it deterministic or stochastic, without wading through the grunt work of finding an analytical solution. It can be used for off-the-cuff estimates or as a proper scientific tool. All one needs to know is how to simulate a given process and its appropriate probability distributions and parameters if that process is stochastic.

Here's how it works:

Construct a simulation that, given input values, returns a value of interest. This could be a pure quantity, like time spent waiting for a bus, or a boolean indicating whether or not a particular event occurs.
Run the simulation a, usually large, number of times, each time with randomly generated input variables. Record its output values.
Compute sample mean and variance of the output values.

In the case of time spent waiting for a bus, the sample mean and variance are estimators of mean and variance for one's wait time. In the boolean case, these represent probability that the given event will occur.

One can think of Monte Carlo Simulation as throwing darts. Say you want to find the area under a curve without integrating. All you must do is draw the curve on a wall and throw darts at it randomly. After you've thrown enough darts, the area under the curve can be approximated using the percentage of darts that end up under the curve times the total area.

This technique is often performed using a spreadsheet, but that can be a bit clunky and may make more complex simulations difficult. I'd like to spend a minute showing how it can be done in Python. Consider the following scenario:

Passengers for a train arrive according to a Poisson process with a mean of 100 per hour. The next train arrives exponentially with a rate of 5 per hour. How many passers will be aboard the train?

We can simulate this using the fact that a Poisson process can be represented as a string of events occurring with exponential inter-arrival times. We use the sim() function below to generate the number of passengers for random instances of the problem. We then compute sample mean and variance for these values.

import random

PASSENGERS = 100.0
TRAINS     =   5.0
ITERATIONS = 10000

def sim():
    passengers = 0.0

    # Determine when the train arrives
    train = random.expovariate(TRAINS)

    # Count the number of passenger arrivals before the train
    now = 0.0
    while True:
        now += random.expovariate(PASSENGERS)
        if now >= train:
            break
        passengers += 1.0

    return passengers

if __name__ == '__main__':        
    output = [sim() for _ in xrange(ITERATIONS)]

    total = sum(output)
    mean = total / len(output)

    sum_sqrs = sum(x*x for x in output)
    variance = (sum_sqrs - total * mean) / (len(output) - 1)

    print 'E[X] = %.02f' % mean
    print 'Var(X) = %.02f' % variance

Running this code yields the intuitive result: we expect 20 passengers on the train. While this is a very simple example, we could easily replace the sim() function with a simulator of any system of arbitrary complexity and get a pretty good estimation of its behavior.

python-zibopt and scip 1.2

2009-09-21T15:54:00.001-04:00

In case you happen to be playing with python-zibopt at all, know that it has been tested and works perfectly well against ZIBOpt 1.2. All you must do is install the new version of ZIBOpt, update your environment variables per the python-zibopt installation page, and recompile.

Special thanks go out to the ZIB folks for maintaining an excellent API across releases. Happy hacking!

python-zibopt initial public release

2009-07-10T18:09:00.001-04:00

I realize it's been a few months since my last post, but the time has been well used. I spent the beginning of this year solidifying my understanding of basic combinatorial optimization techniques in a graduate class taught by Prof. Hoffman at GMU. Being a Linux user and open source advocate, I quickly found that SCIP and ZIMPL were my best options for solving the sometimes rather complex modeling tasks we were assigned. However, we were occasionally given problems that required either going into the C interface for SCIP or doing a significant amount of manual work.

Take, for example, solving a TSP instance using subtour elimination constraints. (See the earlier post about this.) If a TSP is beyond nontrivial size, then the constraints required to describe it can't fit into memory. Thus, one ends up having to solve combinatorial relaxations of the problem repeatedly, adding formulation cuts after each solve until there aren't any disjoint subtours.

With a background in dynamic languages, I wasn't about to dive into a C API just to do this. It began to dawn on me just how much better life would be with a Python interface to SCIP. Among other things, this would give:

An easier time writing code to deal with special problem structures
Access to Python's rich set of libraries for things like CSV and JSON
Proper interfaces to relational databases
Ties to advanced web programming tools for decision support
Multiprocessing (for the truly insane)

Thus, python-zibopt was born. I began coding it up in the days before attending the MIP 2009 Workshop in June, where I had the pleasure of meeting some of the ZIB folks. Since then, python-zibopt has come along fairly quickly. There's a lot to be done, but I've made an initial public release available for anyone who'd like to start playing with it. Here's what it has so far:

Interfaces to SCIP solvers, variables, constraints, and solutions
The ability to turn on and off solver verbosity
Solutions that are false in case of infeasibilty or unboundedness
Mapping from SCIP errors to Python exception types
Solvers can be restarted with additional constraints and variables
Primal solutions can be fed to the solver and are checked for feasibility

Take a look at the examples directory to get started, particularly at tsp.py. This is a TSP solver that can be run against any JSON input file containing a lower triangular distance matrix. It works using the method outlined above. I've included a JSON file with the data from the famous 1954 TSP paper by Dantzig, et al. (If you haven't read this, you really should!)

on the beauty of power sets

2009-02-27T03:52:00.039-05:00

One of the difficulties we encounter in solving the Traveling Salesman Problem (TSP) is that, for even a small numer of cities, a complete description of the problem requires a factorial number of constraints. This is apparent in the standard formulation used to teach the TSP to OR students. Consider a set of $n$ cities with the distance from city $i$ to city $j$ denoted $d_{ij}$. We attempt to minimize the total distance of a tour entering and leaving each city exactly once. $x_{ij} = 1$ if the edge from city $i$ to city $j$ is included in the tour, $0$ otherwise:

$\min$	$z = \sum_i \sum_{j\ne i} d_{ij} x_{ij}$
s.t.	$\sum_{j\ne i} x_{ij}$	$=$	$1$	$\forall i$	leave each city once
	$\sum_{i\ne j} x_{ij}$	$=$	$1$	$\forall j$	enter each city once
	$x_{ij}$	$\in$	${0,1}$	$\forall i, j$

This appears like a reasonable formulation until we solve it and see that our solution contains disconnected subtours. Suppose we have four cities, labeled $A$ through $D$. Connecting $A$ to $B$, $B$ to $A$, $C$ to $D$ and $D$ to $C$ provides a feasible solution to our formulation, but does not constitute a cycle. Here is a more concrete example of two disconnected subtours {(1,5),(5,1)} and {(2,3),(3,4),(4,2)} over five cities:

ampl: display x;
x [*,*]
:   1   2   3   4   5    :=
1   0   0   0   0   1
2   0   0   1   0   0
3   0   0   0   1   0
4   0   1   0   0   0
5   1   0   0   0   0
;

Realizing we just solved the Assignment Problem, we now add subtour elimination constraints. These require that any proper, non-null subset of our $n$ cities is connected by at most $n-1$ active edges:

$\sum_{i \in S} \sum_{j in S} x_{ij} \leq |S|-1 \forall S \subset {1, ..., n}, S \ne O$

Indexing subtour elimination constraints over a power set of the cities completes the formulation. However, this requires an additional $\sum_{k=2}^{n-1} \left(\matrix{n \cr k}\right)$ rows tacked on the end of our matrix and is clearly infeasible for large $n$. The most current computers can handle using this approach is around 19 cities. It remains an instructive tool for understanding the combinatorial explosion that occurs in problems like TSP and is worth translating into a modeling language. So how does one model it on a computer?

Unfortunately, AMPL, the gold standard in mathematical modeling languages, is unable to index over sets. Creating a power set in AMPL requires going through a few contortions. The following code demonstrates power and index sets over four cities:

set cities := 1 .. 4 ordered;

param n := card(cities);
set indices := 0 .. (2^n - 1);
set power {i in indices} := {c in cities: (i div 2^(ord(c) - 1)) mod 2 = 1};

display cities;
display n;
display indices;
display power;

This yields the following output:

set cities := 1 2 3 4;

n = 4

set indices := 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15;

set power[0] := ; # empty
set power[1] := 1;
set power[2] := 2;
set power[3] := 1 2;
set power[4] := 3;
set power[5] := 1 3;
set power[6] := 2 3;
set power[7] := 1 2 3;
set power[8] := 4;
set power[9] := 1 4;
set power[10] := 2 4;
set power[11] := 1 2 4;
set power[12] := 3 4;
set power[13] := 1 3 4;
set power[14] := 2 3 4;
set power[15] := 1 2 3 4;

Note how the index set contains an index for each row in our power set. We can now generate the subtour elimination constraints:

var x {cities cross cities} binary;
s.t. subtours {i in indices: card(power[i]) > 1 and card(power[i]) < card(cities)}:
sum {(c,k) in power[i] cross power[i]: k != c} x[c,k] <= card(power[i]) - 1;

expand subtours;

subject to subtours[3]:  x[1,2] + x[2,1] <= 1;
subject to subtours[5]:  x[1,3] + x[3,1] <= 1;
subject to subtours[6]:  x[2,3] + x[3,2] <= 1;
subject to subtours[7]:  x[1,2] + x[1,3] + x[2,1] + x[2,3] + x[3,1] + x[3,2] <= 2;
subject to subtours[9]:  x[1,4] + x[4,1] <= 1;
subject to subtours[10]: x[2,4] + x[4,2] <= 1;
subject to subtours[11]: x[1,2] + x[1,4] + x[2,1] + x[2,4] + x[4,1] + x[4,2] <= 2;
subject to subtours[12]: x[3,4] + x[4,3] <= 1;
subject to subtours[13]: x[1,3] + x[1,4] + x[3,1] + x[3,4] + x[4,1] + x[4,3] <= 2;
subject to subtours[14]: x[2,3] + x[2,4] + x[3,2] + x[3,4] + x[4,2] + x[4,3] <= 2;

While this does work, the code for generating the power set looks like voodoo. Understanding it required piece-by-piece decomposition, an exercise I suggest you go through yourself if you have a copy of AMPL and 15 minutes to spare:

set foo {c in cities} := {ord(c)};
set bar {c in cities} := {2^(ord(c) - 1)};
set baz {i in indices} := {c in cities: i div 2^(ord(c) - 1)};
set qux {i in indices} := {c in cities: (i div 2^(ord(c) - 1)) mod 2 = 1};

display foo;
display bar;
display baz;
display qux;

This may be an instance where open source leads commercial software. The good folks who produce the ZIB Optimization Suite provide an AMPL-like language called ZIMPL with a few additional useful features. One of these is power sets. Compared to the code above, doesn't this look refreshing?

set cities := {1 to 4};

set power[] := powerset(cities);
set indices := indexset(power);

uncapacitated lot-sizing formulation

2009-02-20T03:30:00.089-05:00

Uncapacitated Lot-Sizing (ULS) is a classic OR problem that seeks to minimize the cost of satisfying known demand for a product over time. Demand is subject to varying costs for production, set-up, and storage of the product. Technically, it is a mixed binary integer linear program -- the key point separating it from the world of linear optimization being that production cannot occur during any period without paying that period's fixed costs for set-up. Thus it has linear nonnegative variables for production and storage amounts during each period, and a binary variable for each period that determines whether or not production can actually occur.

For $n$ periods with per-period fixed set-up cost $f_t$, unit production cost $p_t$, unit storage cost $h_t$,and demand $d_t$, we define decision variables related to production and storage quantities:

$x_t =$ units produced in period $t$
$s_t =$ stock at the end of period $t$
$y_t = 1$ if production occurs in period $t$, $0$ otherwise

One can minimize overall cost for satisfying all demand on time using the following model per Wolsey (1998), defined slightly differently here:

$\min$	$z = \sum_t{p_t x_t} + \sum_t{h_t s_t} + \sum_t{f_t y_t}$
s.t.	$s_1$	$=$	$d_1 + s_1$
	$s_{t-1} + x_t$	$=$	$d_t + s_t$	$\forall t > 1$
	$x_t$	$\leq$	$M y_t$	$\forall t$
	$s_t, x_t$	$\geq$	$0$	$\forall t$
	$y_t$	$\in$	$\{0,1\}$	$\forall t$

According to Wolsey, page 11, given that $s_t = \sum_{i=1}^t (x_i - d_i)$ and defining new constants $K = \sum_{t=1}^n h_t(\sum_{i=1}^t d_i)$ and $c_t = p_t + \sum_{i=t}^n h_i$, the objective function can be rewritten as $z = \sum_t c_t x_t + \sum _t f_t y_t - K$. The book lacks a proof of this and it seems a bit non-obvious, so I attempt an explanation in somewhat painstaking detail here.

Proof:
$\sum_t p_t x_t + \sum_t h_t s_t + \sum_t f_t y_t$	$=$	$\sum_t c_t x_t + \sum _t f_t y_t - K$
1. Eliminate $\sum_t f_t y_t$:
$\sum_t p_t x_t + \sum_t h_t s_t$	$=$	$\sum_t c_t x_t - K$
2. Expand $K$ and $c_t$:
$\sum_t p_t x_t + \sum_t h_t s_t$	$=$	$\sum_t (p_t + \sum_{i=t}^n h_i) x_t - \sum_t h_t (\sum_{i=1}^t d_i)$
3. Eliminate $\sum_t p_t x_t$:
$\sum_t h_t s_t$	$=$	$\sum_t x_t (\sum_{i=t}^n h_i) - \sum_t h_t (\sum_{i=1}^t d_i)$
4. Expand $s_t$:
$\sum_t h_t (\sum_{i=1}^t x_i) - \sum_t h_t (\sum_{i=1}^t d_i)$	$=$	$\sum_t x_t (\sum_{i=t}^n h_i) - \sum_t h_t (\sum_{i=1}^t d_i)$
5. Eliminate $\sum_t h_t (\sum_{i=1}^t d_i)$:
$\sum_t h_t (\sum_{i=1}^t x_i)$	$=$	$\sum_t x_t (\sum_{i=t}^n h_i)$

The result from step 5 becomes obvious upon expanding its left and right-hand terms:

$h_1 x_1 + h_2 (x_1 + x_2) + \cdots + h_n (x_1 + \cdots + x_n) = $ $x_1 (h_1 + \cdots + h_n) + x2 (h_2 + \cdots + h_n) + \cdots + x_n h_n$.

In matrix notation with $h$ and $x$ as column vectors in $\bf R^n$ and $L$ and $U$ being $n x n$ lower and upper triangular identity matrices, respectively, this can be written as:

$\left(h_1 \cdots h_n\right)\left(\matrix{1 \cdots 0 \cr \vdots\ddots \vdots \cr 1 \cdots 1}\right)\left(\matrix{x_1 \cr \vdots \cr x_n}\right) = \left(x_1 \cdots x_n \right)\left(\matrix{1 \cdots 1 \cr \vdots \ddots \vdots \cr 0 \cdots 1}\right)\left(\matrix{h_1 \cr \vdots \cr h_n}\right)$

or $h^T L x = x^T U h$.