**Contact**:
iaindunning at gmail,
@iaindunning,
github/IainNZ

MIT has a wonderful thing called IAP - Independent Activities Period - over the gap between Christmas/New Years and the resumption of normal classes in early February. This year I was involved in a couple of things that were a lot of fun.

While there is no question that everyone in the Operations Research Center here at MIT is very talented in the area of mathematics (and many other things!), there is a wide variety in the level of ability when it comes to programming. The main reason for this is probably the variety of paths people take to reach the ORC, from pure mathematics to engineering. This all applies to the faculty too!

Being problem solvers, we decided that one way we can address this issue is for students to transfer their skills to their fellow students. Its hard to fit hands-on, programming-intensive classes into the traditional semester-long lecture framework, so we made our own course that runs in IAP - with this year’s class being the first.

The format was a series of 3-hour-long workshops. The topics covered were

- introductory and advanced classes on R (Allison O’Hair, Andre Calmon, John Silberholz),
- data visualization (Virot Chiraphadhanakul),
- customizing integer programming solvers (Ross Anderson),
- convex optimization (Vishal Gupta),
- distributed computing for optimization (John Silberholz),
- and using Python for mathematical modelling (which was my workshop.)

I really enjoyed delivering my workshop. It was the first time I had given such a long presentation, and it was hard to estimate how long the material I created would take to work with. My guess was we might actually finish short of the 3 hours. My general approach to the workshop was to assume that a) the class was very familiar with linear programming, and b) they knew a programming language, not necessarily Python. I began with a brief introduction to what a modelling language is and the options out there. This was followed by roughly 30 to 45 minutes on Python-specific features - lists, tuples, dictionaries, and list comprehensions. I followed this up with an exercise where I asked the following three questions:

- Write a list comprehension to generate the square of the numbers between 1 and 10
- Write a list comprehension to generate the odd numbers between 1 and 20
- Write code to solve the FizzBuzz question.

Students worked on this by themselves for a few minutes as I walked around the room. Once I sensed that a majority had finished the first two, and some had solved FizzBuzz, I coded the solutions live on the projector.

Next I introduced PuLP itself, the modelling language we used. I demonstrated a couple of common patterns and showed what a simple knapsack problem would look like in PuLP. This was followed by another exercise where I described the “p-median”/facility location problem, provided the mathematical formulation, and asked the class to implement it. The reason I selected this problem is that it is very light on data - I wrote a line to generate the customer locations, and provided that as a template, but otherwise gave them a blank slate. People rose to the challenge, and did a great job on this exercise. In a way this reflects the intuitiveness of PuLP/Python, especially if you have a deep understanding linear programming. Once again I finished this exercise by coding my solution live on the screen.

The next section was focused on how to bring data into the mix, including random number generation and reading from CSV files. This was going deeper into Python, and sometimes introduced a lot of concepts quickly (file IO especially.) I ended up tabbing-out of the slides and going into Python to make it more concrete. The final exercise was basically the classic diet problem, with the data being stored in CSV files. There were many possible ways to store the data in memory, and it was a treat to see the diversity of approaches.

I learnt an interesting lesson here: many chose to store the tabular data in a “dictionary of dictionaries”, which is certainly an intuitive way to do it. Unfortunately the idea of nesting dictionaries was a bit of an escalation for some people, and many of the less confident students felt obliged to do it this way when they heard others talking about it. Although I think most people got it in the end, I probably needed to provide a few different alternatives before hand using different data, to give the students some more guidance. On the other hand, I like that people were able to come up with approaches by themselves. I suppose this is always a tension in teaching.

I finished with some performance concerns and pointed students towards solver-specific interfaces if they need them. There was one final cautionary example that I used and I want to reproduce here. Consider this situation:

- We have lists
`SERVICES`

and`LOCATIONS`

- We have a list
`PROCESSES`

, and for each process, we have the service it belongs to. - We are deciding the assignment of
`PROCESSES`

to`MACHINES`

,`x[p,m]`

- Imagine we have some constraint like the following:

```
for s in SERVICES:
for l in LOCATIONS:
prob+= ... <= lpSum(x[(p,m)] for p in PROCS if PROC_SERV[p] == s for m in MACHS if MACH_LOCATION[m] == l )
```

- Can you spot the issue? The problem is we check every possible process-machine pair for every constraint. Say there ~N services, ~N locations, ~N processes and ~N machines. We go around the loops N^2 times, and the list comprehension does N^2 checks - so N^4 work!
- The “right way” to do it is to calculate a mapping of services to processes once, and a mapping of locations to machines once, e.g
`SERV_PROCS[s]`

and`LOCATION_MACH[l]`

. This will take N^2 work, but you only do it once.

```
for s in SERVICES:
for l in LOCATIONS:
prob+= ... <= lpSum(x[(p,m)] for p in SERV_PROC[s] for m in LOCATION_MACH[l])
```

- Now we do far less than N^2 work per constraint, in fact we do the minimum possible. This makes a huge difference in practice - be on the look out for it!

We ran out of time right after that, so we didn’t have to use my final exercise which was far more involved and could soak up a large amount of time. In fact, I had to rush this last example a bit, so next time I’ll need to go just a bit faster on the exercises, or trim elsewhere. I think the class was a success though, and the feedback was really positive, which was a great feeling - teaching might be addictive!

I’ve mentioned Julia on this blog previously - it is a up-and-coming programming language with a lot of promise. It is expressive, fast, good for math/science, and has some advanced tricks that make it a good candidate to create domain-specific languages. There is already some great package support, and many tasks one might do in R or MATLAB can be done in Julia already.

My partner-in-crime Miles Lubin and myself decided to tackle the challenge of bringing Operations Research tools to Julia. Miles has done some great work on making an interface to the COIN CLP solver and we’ve worked together to make a modelling language that we call JuMP. JuMP stands for “Julia for Mathematical Programming” and it can be compared with PuLP and Pyomo. While making an easy-to-use modelling tool was of course a goal, what we really wanted was to be as fast as commercial tools like AMPL, or maybe even a C++ solver interface. A common situation for OR people to find themselves in is that they test algorithm ideas in a language like Python or MATLAB, but find that it doesn’t scale up for “real problems” of a larger size, forcing people to change to AMPL or C++ to get the performance they need. This is not efficient - our hypothesis was that Julia can both be useful for experimental purposes and production purposes.

We have largely achieved this at this point, and we are on the same speed level as AMPL (within a factor of 2) and far faster than Python-based languages. The real challenge now is to flesh out the module while maintaining the speed, and to build up a body of examples. Other work is required to hook Julia up to other solvers, but this is a somewhat separate issue.

We presented our work on this, and the work Miles did with an efficient simplex algorithm implementation, at a IAP workshop on Julia. It was very interesting to see the work being done, and I have a lot of confidence in its future. Check it out today!

© Iain Dunning, 2015. Base theme/CSS by Skeleton.