Biarri Blog

An efficient and effective research environment (on linux, with vim, jabref and more)

2012-01-25T12:01:00.000+11:00

So, I would like to share the environment that I have created for the purposes of doing research. Specifically it is an environment that allows me to:

Gather research papers,
Comment on them in various ways, and review these comments at large,
Store this information in source control for the purposes of sharing between my machines, and
Write up and deal with ideas in a systematic fashion.

So, the perhaps the first component of this system is, what operating system? Happily, it doesn't exactly matter. I use Ubuntu, but I also use Windows 7. The great thing about this scheme is that it is adaptable to any environment that runs the tools (well, obviously) and the tools all have multi-platform support.

I personally find editing in vim nicer on Ubuntu, and there is one or two arguably minor things that linux has that Windows does not (XMonad, for example), but I will elaborate on these later.

The general scheme

This approach is, obviously, tailored specifically for me, and given that I have a significant programming background, I am happy to solve some problems with actual programming. I also quite enjoy "systematic" approaches; i.e. things that may not neccessarily be the most user friendly or easy, but ones that follow a specific and clear pattern that makes logical sense.

This approach may not suit everyone, but hopefully there are at least interesting and useful ideas here that you could adapt to your own situation.

The beginning - reference management

So, of course in order to gather research papers it is neccessary to store them in a useful way. JabRef is free, and is a very nice option for this. I've described my custom JabRef set up on my own personal blog a few months ago, please read that for how to do this.

One thing on using JabRef is that sometimes you need to correct the format that bibtex exports give you. For example, one thing I am often changing is lists of authors like: "Billy Smith, Joey Anderson" to "Billy Smith and Joey Anderson"

It's not immediately clear to me why the bibtex are generated the wrong way from some of these sites, but nevertheless. This simple correction is neccessary for the data to be stored properly, and the authors to be picked up correctly, etc.

Okay, now that you've read that, you understand that I save all my PDFs to a specific location, a folder like ".../research/reference library". Where is this folder? It's in the research repository.

The "research" repository

I keep all my research, on any topic, in one generic folder, called "research". This is a private git repository hosted on BitBucket.org. I chose bitbucket over github because bitbucket has free unlimited-space private repositories, while githubs cost money. It is neccessary for the research repository to be private for two reasons, one obvious one is that it contains paywall-restricted PDFs, and the other is that it's just not appropriate to have in-progress research notes be viewable by anyone, mainly because people may see how silly my comments really are!

So, the general structure of my research repository is as follows:

~/research
    /articles
    /conferences
    /diary
    /jabref
    /reference library
    /projects
        /semigroups
        /...
    /quantum-lunch
    /...

The contents of these folders are as follows:

~/research/articles - Articles

This contains folders that map directly to papers that I'm trying to write (or more correctly at the moment, scholarships that I'm applying for, and misc notes). These are all, unsurprisingly, LaTeX documents that I edit with Vim.

When I complete an article, I created a "submitted" folder under the specific article, and put the generated PDF in there, but up until that time I only add the non-generated files to source control (this is the generally correct practice for using any kind of source control; only "source" files are included, anything that is generated from the source files is not).

~/research/conferences - Notes on conferences

In here I have folders that map to the short conference name, for example "ACC" for Australian Control Conference. Under that, I have the year, and within that I have my notes, and also any agenda's that I may have needed, to see what lectures I would attend. The notes should be in vimwiki format (I will describe this later) for easy/nice reading and editing.

~/research/dairy - Research diary and general ideas area

This is the main place I work on day-to-day. It contains all my notes, and minutes from various meetings and lectures I attend. It contains a somewhat-daily research diary, and a list of current research ideas, past research ideas (that were bad, and reasons why) and so on.

My preferred note taking form is vimwiki (to be described below), so in here are purely vimwiki files.

It's not essential that you also use vim (and hence vimwiki), but it is appropriate that whatever mechanism you use, it is a format that is ameneable to source control (i.e. allows nice text-based diffs). Emacs or any plain-text editor will be sufficient here.

~/research/jabref - Bibtex files

This is perhaps not the most appropriately named folder, but nevertheless. It contains all my .bib databases. I actually only have 3. One is very inappropriately called "2010.bib", with the view that I would store research by the year I gathered it. I'm not following this approach and I actually just keep all my research related to quantum computing (and more general subjects) in here.

I have two other bib files, one is related to a secondary field of that I am interested in researching. That is to suggest, in 2010.bib I have only documents related to quantum computing, theoretical physics and some theoretical computer science. I have a different .bib for research in completely seperate fields, say investment. The other is "lectures.bib", and it is obvious what that contains.

It's worth noting that I actually have two systems for storing lectures. One is the above, where the lecture set fits into a nice single PDF. The other is when the lecture series is split over several PDFs. These I store under a generic "University" repository that I use for my current studies (including assignments and so on). This component of my current setup needs work, and I'm open to suggestions here.

~/research/reference library - All the PDFs

Every PDF releated to my research is stored in here, prefixed with the bibtex key. So, for example I have "Datta2005, 0505213v1.pdf". JabRef associates this with the relevant item in the .bib file by the prefix, and virtue of this link in the .bib I have a trivial way to programmatically (if I'm so inclined) get this PDF.

I don't ever browse this folder directly, and currently it contains ~600 PDFs and is about 1 gig in size. Storing this much data in a git repository may offend some people, but essentially they are wrong to be offended. It is okay to store binary files as long as they are not constantly changing, and they are considered a key component of the repository; which in this case they are.

There are some downsides to this, though, and I think it's plausible to consider alternative arrangements.

The viable alternatives are:

Dropbox for PDFs,
Secondary git repository for PDFs only, and include as submodule, or
Some other non-local storage (say, the one offered by Zotero).

I dislike all of them because for me I prefer to have everything together. I could see dropbox being suitable, because technically it's not neccessary to have verioned PDFs.

If you have any comments on this, please share them.

~/research/projects - Specific Projects

You'll notice I have one folder here, "semigroups". This relates directly to a research scholarship I completed at RMIT. This actually involved a python program, which I have in this directory, as well as some miscellaneous files. It may be appropriate to have nicer codenames for projects, or somehow related them directly to the scholarship details. I think the best approach here is to have a codename, which is detailed in the "diary" folder, and then there is no risk on confusion or duplicate names. The scholarship details could be held seperately in the folder, because perhaps the work could be continued across scholarships.

Anyway, it's probably not neccessary to overwork this structure. It can always be changed, and it shouldn't be prohibitively difficult.

~/research/quantum-lunch - Files related to my reading group

This folder is indicative of the other types of folders you may find in this directory. In here, I have some misc python scripts related to this group. There are no notes in here, they are kept in the "dairy" folder.

Technically this should be a transition area, where scripts and programs that reach an appropriate level of maturity/usefulness are either published publically (in a different repository), or moved to an appropriate folder under project, but I've not yet gotten to that stage.

It's worth noting that I do have a public github profile: silky, underwhich I will, and do, publish and tools that are worthwhile being public. If one of these projects reaches that stage, I'd essentially move it out of here (delete it) and continue it in that area.

The tools

So, with the repository layout described, let me know discuss the tools I use. We've already covered JabRef, for reference management, so we have:

JabRef (as mentioned), for reference management,
Vim + Vimwiki plugin, for taking notes, keeping ideas, and writing LaTeX,
Okular, for reading PDFs, and annotating them [linux],
Python, for programming small scripts,
XMonad, for window management [linux], and
pdflatex and bibtex, for compiling latex (from the TeXLive distribution).

So, almost all of these are available on any platform. Okular is worth commenting on, because it has an important feature that I make use of - it stores annotations not in the PDF but in a seperate file, that can then be programmatically investigated. If you can't use okular, then you may find that your annotations to PDFs are written back into the PDF itself, and it will be difficult to extract this. You can decide whether or not this bothers you when I describe how I use my annotations.

I will now describe the usage pattern for the various tools, starting in order for easiest to hardest.

Tools - Okular

So, install okular via your favourite method, say "sudo apt-get install okoular", and then open it. You will want to make it your default PDF editor, and I also choose to have it be very minimal in it's display; setting the toolbar to text only, hiding the menu, and hiding the list of pages on the left. I also configured a shortcut for exiting, namely pressing "qq".

For me this is indicative of an important principle - make small customisations that improve your life. It's worth thinking about, as they can often be trivial, but provide a nice noticable benefit.

You will also want to enable the 'Review' toolbar. This allows you to highlight lines of interest, and also add comments. Your comments are saved in a location like:

~/.kde/share/apps/okular/docdata/[number].[pdf-name].xml

This is where it gets fun. I've written a program to capture these comments, as well as comments in the 'Review' field of the .bib file. This tool is available on my github: get-notes.

You may need to adjust the 'main.conf' to suit your needs, or even change the source in some fashion. The code is pretty trivial, but requires some python libraries that you can install with easy_install.

This file products vimwiki output (you can trivially change this however you like, if you program in python). I then symbolically link this generated file ("AutoGeneratedNotes.wiki") to my "~/research/diary". Of course, following the general strategy of not including generated files in the source code, I do not break this rule for this file. There is one perhaps obvious downside to this: The output might be different on different machines, because the ~/.kde/... folder is not under source control. I consider this acceptable, because this file is a "transitional" file, in that it is not supposed to be for long-term storage of ideas and comments.

The contents of this file should be reviewed, occasionally, and then moved into either a PDF of comments, or into the research diary for an idea to investiage, or removed because you've dealt with it.

For example, I have a comment in the "Review" field of the file "Arora2002". It says: "Has a section on the Fourier Transform that might be interesting". This should, eventually, be transitioned into a minor topic to investigate further, or a small writeup in a private "Comments on things" LaTeX document, where you write up, slightly more formally, and with maths, your thoughts on things you've learned. I have this document under my "articles" folder.

With this in mind, it is then not an issue that the generated output is different bewteen machines, because ideally there will be no output on any machine, one it has been sufficiently transitioned.

Tools - LaTeX

As indicated, I use LaTeX to write up maths and more detailed notes, proposals, applications, etc. You may wish to use some front-end for LaTeX authoring, for example LyX, but as I already do a lot of work in vim, I prefer to also do LaTeX in here. If I were to switch to another editor, it would probably be Emacs.

Tools - Python

As mentioned in the above comment, I use python to write small scripts. Because they're in python, they are essentially directly runnable on any system (provided the associated packages can be installed).

I also like python because it provides various packages that are actually useful for my research (like, for example, numpy). You can get similar funcionality from free Math environments, though, such as Octave.

Tools - XMonad

XMonad is not particularly neccessary for this workflow, but I include it because I find it's ease of use aids in efficient reading and editing. I don't want to go into significant detail of XMonad configuration (but it's a fun way to spend your time), you may simply review my XMonad configuration on github.

What I like about it is the concept of focus. You can simply and easily make a PDF full screen, for distraction-free reading, and then switch things around to have vim side-by-side for commenting with context.

Feel free to disregard this, if you are using Windows, as it's equally possible to do fullscreen and side-by-side editing in Windows 7. XMonad also offers other benefits for general programming, which is the main reason I have it.

Tools - Vim + Vimwiki + LaTeX Box

Essentially the last item in my setup is Vim. It's hard to express the level of obsession one has for Vim, after a while of using it. It is highly customisable, and includes and inbuilt help system, which I used all the time, when initially learning it.

Most people will find Vim initially difficult to use (I did, when I first learned it when starting work here), but if you dedicate a few days to using it correctly, and you make significant use of the ":help [topic]" command, you will get the hang of it.

You aren't truly using Vim correctly (or, indeed, living a full life), if you don't get various plugins. The neccessary ones for LaTeX + note taking are: Vimwiki and LaTeX Box or Vim-LaTeX.

You can find the current state of my Vim configuration, again on my github - .vim

I actually currently use Vim-LaTeX, but I am planning on changing to LaTeX Box because it is more lightweight, so I would recommend starting with LaTeX Box.

The nice thing about using okular is that you can recompile your LaTeX document with the PDF open, and it will refresh it, keeping the current page open. This is very useful when typing long formulas, and reviewing your work.

I have configured Vimwiki to open with "rw", so I can type this at any time, in Vim, and be looking at the index to all my research notes. In this I have links to all my diaries, my storage spots for old search ideas, and a big list of topics to look into. I also make "TODO" notes in here, and review them with one of my other tools, "find-todo" (on the aforemention github, under /utils). This gives me a list inside Vim, and I can easily navigate to the appropriate file. Again, the TODO's are items that should be transitioned.

Review

I have documented my reseach environment, as it stands currently. It allows me to make notes easily, transition them in an appropriate workflow, and access all my documents at any time, from any computer.

The proof of a good research environment obviously in the blogging, it's in the producing of legitimately good research output, and of course that's yet to be delivered (by myself personally), so it's not possible to objectively rate this strategy for it's actual effectiveness. Nevertheless, I do feel comfortable with this layout; I feel like I can take the appropriate amount of notes; I feel my notes are always tracked, and I feel that I have a nice and readable history of what I've done. I like that I can track bad ideas; I like that I can make comments "anywhere" (i.e. in Okular or in JabRef) and have them captured automatically for later review, and I like the feeling of having everything organised.

I hope this description has been useful, and I would love to hear about any adjustments you'd propose, or just your own research strategies.

--Noon

Vim (and Emacs)

2011-12-15T16:57:00.001+11:00

Biarri now has a few dedicated vim users with their configuration files available online: Myself (Loki) and Noon.

My configuration file is based on Sontek's with a few changes and additions. His description at http://sontek.net/turning-vim-into-a-modern-python-ide is worth a read.

You'll notice that I have disabled my arrow keys. Despite liking the Cure and The Smiths this isn't to hurt myself. It's not ever so I keep my hands on the home row and use hjkl. Instead it's to force me out of insert mode. I now stay in normal mode as much as possible. Since disabling my arrows I've found the combination of traditional vim motions (especially f/F) combined with Easy Motion is fantastic. I usually use f for short jumps and Easy Motion for longer movements.

Easy Motion is really best explained by watching the demo video closely, or just using it. It's quite wonderful and works well with visual mode. I'm also enjoying Command-T, PyFlakes, snipmate and ack. I'm still confused about what the best autocomplete solution is. Supertab doesn't seem to work right for me and ACP, which I was using for a while is useful but somewhat annoying. I would like it to show call signatures and if requested, docstrings in a non-annoying way too. Command-T is the greatest thing imaginable for opening files, buffers and tags.

I've also installed vimwiki but i'm so far unsure about it. I used it to store my ideas and todos in a more organised manner. I suspect it's a bit like a less awesome org mode. Speaking of which we now also have an Emacs user: https://bitbucket.org/tobyodavies/shared/src If anyone knows a good python autocomplete solution I'd love to hear about it.

Loki

Subscription payments

2011-12-09T17:36:00.001+11:00

Selling things on the internet is quite easy and there are plenty of way to process credit cards and accept online payments. However subscriptions / recurring payments aren't easy. They aren't easy at all. Subscription changes, card failures, multiple subscriptions for a single user (different workbench products in our case) and multiple currency billing make it hard. Dealing with banks 80's style systems makes it even harder.

When we were first working out how to cope with online payments we found lots of good looking options, all of which were only available to American companies or via paypal. (Recurly, CheddarGetter Spreedly) Eventually we found RBS Worldpay.Worldpay turned out to be too inflexible and far to much effort to set up. The time required was also insane. From the start of our application to ready to get our payments up took 7 months and significant amounts of cash. We gave up on Worldpay when we were trying (hard) to integrate it. Their system for recurring payments is so ill documented and inflexible to be useless. Users appear to require a separate login and password for Worldpay. This is unacceptable. Currently we're looking into Saasy, from fastspring. They respond fantastically quickly to support request and got our testing site up and approved us in a day. So far it looks great, flexible and easy but I'll post another blog when we've got it all live. I'll also put some code for the server side part in python using cherrypy up on github, when it's a bit more polished.

Loki

Gartner’s top-10 strategic technologies for 2012

2011-12-01T08:22:00.000+11:00

Gartner’s top-10 strategic technologies for 2012

Selecting a cell in an HTML table with jQuery

2011-11-01T14:46:00.005+11:00

I wanted to get (and set) a cell value in a table using jQuery, selecting it by row and column index (its x and y position in the table, if you like). Frustratingly, I couldn't find a simple one-liner to do it: I could find code that gets all the cells, or get them by ID (but I don't want to have to make explicit IDs for each cell of my table), or get all the values in a given column (close, and did lead me in the right direction).

So here's two similar ways to select a given cell, where row and col are my x, y positions, and my_table is the ID of the table:

$('#my_table tr:eq('+row+') td').eq(col)

and

$('#my_table tbody tr:eq('+row+') td:eq('+col+')')

PostgreSQL Tricks

2011-10-18T15:29:00.003+11:00

Over the last year or so I've gathered a little trick bag of PostgreSQL recipes. Here are some of the best.

These all work in PostgreSQL version 8.4.7 and probably most other versions. The "--" below are PostgreSQL comments.

To see your PostgreSQL version, use:

select version();

It's sometimes useful to be able to discover the tables in a schema automatically. To do this you can use the following command:

>>> \dt my_schema.* -- here the * is a wildcard

Or this command, which uses the PostgreSQL internal tables pg_class and pg_namespace:

>>> select n.nspname, c.relname from pg_class c join pg_namespace n on n.oid=c.relnamespace where n.nspname='my_schema';

Where "my_schema" is a schema name. Once you've found the table you're interested in, \d table_name gives you the columns and their types.

To discover the columns in a table (that is in a schema):

>>> select a.attname from pg_class as c, pg_namespace as n, pg_attribute as a where n.oid=c.relnamespace and n.nspname='my_schema' and c.relname='my_table' and a.attrelid=c.oid;

Where "my_table" is the table name in your schema. This uses the PostgreSQL internal tables pg_class, pg_namespace and pg_attribute.

To discover the type of a column in a table (that is in a schema)

>>> select a.attname, t.typname, a.atttypid from pg_class as c, pg_namespace as n, pg_attribute as a, pg_type as t where n.oid=c.relnamespace and n.nspname='my_schema' and c.relname='my_table' and a.attrelid=c.oid and a.attname='my_column' and t.oid=a.atttypid;

Where "my_column" is the column name in the table. This is similar to the previous command, but also uses internal table pg_type.

If you have some duplicate rows in your table, you can delete them using the hidden ctid column, e.g.:

>>> select ctid, * from my_table; -- show me the repeats

>>> delete from my_table where CAST("ctid" as text) like '%3%'; -- kill some according to some wildcard pattern

A more automatic way is:

>>> delete form my_table where ctid not in (select max(dup.ctid) from my_table as dup group by dup.x, dup.y, dup.z); -- substitute appropriate column names for x, y, z

Allowable formats for dates are controlled by the PostgreSQL setting 'datestyle'. To see it, type:

show datestyle;

To change the datestyle you can use (for example):

set datestyle to 'iso, dmy';

But note that this change is not a permanent change (it only applies to this session/cursor).

Lastly, I sometimes use this command:

table my_table;

as a shorthand for select * from my_table.

Debugging Internet Traffic

2011-10-14T12:12:00.002+11:00

Below are some tools and thoughts that may be useful to people who are having difficulties trying to debug issues related to internet traffic.
There are a number of other tools (e.g. Firebug, Wireshark/ethereal, etc) that are not mentioned below. However I thought a quick one pager on some of the tools that I have used and found useful might be useful to someone else.

Tamper (Firefox plug-in)

https://addons.mozilla.org/en-US/firefox/addon/tamper-data/
This is a Firefox plugin which allows users to view the exact data they're sending out through their browser. It's useful for debugging traffic, as well as those times when your browser has remembered your password and you've forgotten your password (as you can see the password that is being sent). It also gives you the opportunity to intercept this data and change it which can also be useful for debugging purposes. The disadvantage with this tool though is that it will only work on data sent through Firefox. Fiddler2 (below) can be used to intecept traffic from a wider range of sources (e.g. iexplorer, firefox, wget, curl, ... and other applications).

wget and curl

http://www.gnu.org/s/wget/
http://curl.haxx.se
These tools are very similar and are useful for sending web requests via command line. They both work in Linux and Windows (and other OSes). I was using this to help debug an issue we were having at a client site that was not allowing our newt traffic to be sent to our servers. It means you can easily change things like proxy settings, username/passwords, URLs, the posted data - without having to recompile an exe and copy it onto the client machine. I personally found curl worked better than wget. I couldn't get the username/password information to be sent correctly with wget. There was plenty of help on the net, but it just didn't seem to work for me for some reason. A typical example of the types of curl requests I was sending looked like this:
curl -L -x proxy_details:proxy_port --user username:password --output outuptfile.html http://website.com/method --data "user_info=stuffhere%3Bmorestuffhere%3Bextra_info=morestuffhere"
These command line tools can of course called via scripts as well if required.

Filddler 2

http://www.fiddler2.com/fiddler2/
I found this tool paricularly useful when diagnosing firewall/proxy issues as I could see exactly why various requests from curl/wget were, or were not, working.

I'm certainly not an expert on Fiddler2, and I'm not abou to try and tell you all the ins and outs of how to use this tool: firstly because I don't know all of them, secondly because a quick Google search will help you there.
However, some of the useful features that I have found useful are:

It can monitor traffic from differnt sources. The "Process" column in the left pane in the above image shows that curl and firefox traffic has been detected.
The upper pane on the right shows the request details of the selected item in the left pane:
- Headers
- TextView
- WebForms
- HexView
- Auth (useful if want to see how usernames and passwords are being handled)
- Raw
- XML
- JSON
The lower pane on the right gives the response details of the selected item in the left pane:
- Auth
- Caching
- Privacy
- Raw
- XML
- JSON
Another one of the cool features of fiddler2 is that it can be used to debug encrypted traffic. It does this by installig a security certificate on your computer. Ideally you would only install this security certificate on a machine for debuggnig purposes and one that was not passing sensitive information over the internet as it poses a security risk. It should be removed from the computer once you have finished debugging (it's easy enough to put another back in the future if requred).

Kevin Slavin: How algorithms shape our world

2011-10-13T16:08:00.002+11:00

Pre-emptive optimisation; or, how to be deliberately wasteful of server resources

2011-09-28T08:59:00.002+10:00

For one of our long-standing clients we have been running vehicle routing optimisations on a daily basis. A file of daily orders is uploaded into our Workbench system, and is split up into several regions, each of which needs to be separately optimised. A planner works through each region, going through a series of data checks (e.g. location geocode checking), before hitting an “Optimise” button.

All of the heavy lifting is done on a server (the planner accesses it through a web app via a browser), so it’s possible that the server could silently start up the required optimisations without the planner’s involvement, and in the (fairly common) case where the region does not require any data corrections, when the planner is up to the optimisation stage, the result could be immediately available (as it has already been run, or is in the process of being run). This idea now been implemented and took only a short amount of Python code.

Furthermore, it runs in parallel as each optimisation is itself split into a separate child process (running a C++ exe) which Linux distributes across the 8 cores of our server machine.
The pre-emptive optimisations are kicked off using Python’s multiprocessing package, as follows:

from multiprocessing import Process

...

p = Process(target=start_preemptive_optimisations, args=(...))

p.start()

Control is returned to the user at this point while the optimisations run in the background. Results are stored in a folder whose name is stored in a database table; when the planner then comes to press Optimise, the system checks if there have been any data corrections – if so, it runs the optimisation from scratch as usual (the pre-emptive result for that region is thus never referenced); however, if there are no corrections, the system simply uses the stored result from the folder.

The end result for the user is that in many cases the optimisations appear to run almost instantaneously. There are really no downsides as we do not pay for our servers on a CPU cycle basis, so we can easily be wasteful of our server CPU time and run these optimisations even if their results are sometimes not needed.

One “wrinkle” we discovered with this is that we had to make our process checking more robust. There is Javascript in our browser front end that polls for certain events, such as an optimisation finishing, which is indicated by a process ceasing to exist. The Python for this is shown below, where “pid” is a process ID. The function returns True if the given process has finished or not.

def check_pid_and_return_whether_process_has_finished(pid):
if pid and pid > 0:
multiprocessing.active_children() # reap all zombie children first; this also seems to pick up non-children processes
try:
os.waitpid(pid, os.WNOHANG) # this reaps zombies that are child processes, as it gives these processes a chance to output their final return value to the OS.
except OSError as e:
if int(e.errno) <> 10: # the 10 indicates pid is not a child process; in this case we want to do nothing and let os.kill be the function to throw an exception and return True (indicating process is finished).
return True
try:
os.kill(pid, 0) # doesn't actually kill the process, but raises an OSError if the process pid does not exist; this indicates the process is finished. Applies to all processes, not just children.
except OSError as e:
return True
else:
return True
return False

Note the “reaping” of zombie processes here, and the complication that arises if a process is not a direct child of the calling Python process (it might be a “child of a child”). In this case we use a call to the (in this case rather mis-named) function os.kill.

Spatial Clustering in C++ (post 5 of 5) – The final outputs

2011-09-05T10:49:00.010+10:00

In this fifth and final post of the series on Spatial Clustering in C++, we’ll look at some of the final outputs produced and wrap up the discussion of techniques and codes used.

In the last post I presented some code that can be used to generate an Alpha Shape from a set of points, create it as an OGR Polygon, and output these polygons to MID/MIF files. I’ve since found that one of the functions used in that code, GetNextRing, is quite inefficient when the polygon has a complex perimeter (many points on its boundary), which happens with very small values for alpha. Still it wouldn't be too hard to speed that code up.

It’s interesting to generate some different Alpha Shapes with varying of the alpha parameter. Recall that alpha corresponds to how tightly the polygon is “shrink wrapped” around the points it encloses; it is actually the radius of circles that define the “bites” taken out around the perimeter of the shape. If alpha is too high, one gets a very approximate shape that starts to approach the convex hull of the points, while if alpha is too low, the shape starts to get very “spiky”, and eventually one starts to lose whole parts of the polygon completely. The good news is that the polygons produced are not too visually dissimilar if you stay within reasonable bounds for alpha – that is, they are not overly sensitive to its value.

Let’s have a look with the San Francisco bay area as our example. Here I’m using an eps in the DBSCAN of 0.009 with minPts = 8, which gives the following clusters:

And for the Alpha Shapes, an alpha of 2500 gives the following:

In the next picture you can see how the shapes represent the points (note also the hole in the large green cluster – the Alpha Shapes code we’re using doesn’t allow for actual holes within the shapes produced):

Here I’ve used Bing Aerial in MapInfoPro to give the shapes produced some spatial context, and I’ve used 33% transparency on the polygon layer (use Layer Properties). Incidentally you can use Style Overrides to give patterned output as well, such as:

With alpha 5000, the results are very similar, though some of the “holes” along the polygon boundaries have now been “filled in”:

The colours of the shapes are different here just because of the funky unique RGB generator function I’m using which depends on the creation index of the polygon; the clusters are the same, only alpha is different.

Even with alpha = 20000, the shapes are still good visual approximations of the points:

With 50000 the shapes are still OK but some big holes have now been filled in:

And with alpha 5000000 one starts to get nearer to the convex hull:

With smaller values of alpha, such as 1000, the shapes are similar to alpha 2500. Below that, however, the Alpha Shapes algorithm starts to produce multiple, very spiky (many points defining a complex boundary) polygons for each cluster, such as this example with alpha 500:

These anemic shapes no longer approximate the points very well.

With larger data sets just showing the clustered points gets slower, so you could potentially also make a geoset and use zoom levels, with polygons shown at higher zoom levels, then show the clustered points when zoomed in closer.

In terms of output, I’ve found Tab2tab is a useful little program to convert between Mid/Mif and Tab formats, though it is not without its problems: older versions have a loss of precision, and current versions appear to lose the colours of the polygons (regions in the parlance). One could potentially also write out the polygons or clusters using OGR functions directly in the output format desired. And if you don’t have MapInfoPro, QGIS is a good open source alternative.

If you've made it this far, then I hope you've enjoyed this five part blog post. I'll leave you with a summary of the tools and algorithms that have been deployed in this exercise:

Tool/algorithm	What it does	Notes/shortcomings
TIGER spatial data	Source for free US spatial data	In this exercise, used as a data source of points for spatial clustering
MapInfo	Display and edit spatial data	Commercial product, but very good
kdtree	Data structure for storing 2 dimensional points	Allows fast lookup of nearby points. Good C++ implementation here.
OGR	Library for computational geometry	Good utility functions for reading and writing points and polygons
DBSCAN	Algorithm for spatial clustering	Simple to implement; previous blog posts in this series show how to make this run really fast.
Convex Hull	Shape containing a set of points	Not a good approximation to spatial clusters, but easy and fast to compute e.g. with OGR functions
Clarkson code for Alpha Shapes	Old style C code for producing Alpha Shapes (polygons representing sets of points)	Slightly dodgy code, doesn't produce polygons with holes inside, but runs OK if called as a black box executable.
Tab2Tab	Map data conversion utility	Older versions lose precision; new versions lose polygon style/colour information
QGIS (Quantum GIS)	Open source alternative to MapInfoPro	Slower and not as feature rich as MapInfoPro

Spatial Clustering in C++ (post 4 of 5) – Generating Alpha Shapes

2011-08-30T16:30:00.008+10:00

In previous posts we saw how to run the DBSCAN spatial clustering algorithm over a set of spatial data points: in our example, all block centroids in California. We also saw how to make it run really quickly. In this post we’ll investigate a way to turn those clusters into polygons. These polygons will represent our clusters in an approximate way, but will be fast to display and provide a visually compelling way to appraise the clusters we’ve created.

One way to get a polygon from a set of points is to form the Convex Hull. However, for our clusters, convex hulls, though easily generated using OGR’s ConvexHull() on a MultiPoint object, do not look so great – they cover over large swathes of what we’d intuitively think of as the “shape” describing the points. An example from a real cluster produced by DBSCAN is shown below:

In this example, the light green coloured cluster would be badly represented by a convex hull (which is shown in pink):

Another example can be found in Part 1 of this blog series, for the largest cluster in the picture.

Unfortunately there’s no unambiguous way to wrap a set of points in a non-convex polygon (that is, there’s no such thing as a “bounding polygon” – think about it...). However, we can use something called Alpha Shapes, which have a parameter (alpha) which describes “how concave” we’d like our shape to go (intuitively, how tightly to “shrink wrap” our points). Good descriptions of Alpha Shapes are provided here, here and here, and in the context of spatial data, here.

For the above example, with a suitably chosen alpha (in this case, 2500), the alpha shape looks like the following:

There’s not a lot of Alpha Shape C++ code around and it’s not something you’d probably want to code up yourself unless you’re very interested in computational geometry and have lots of spare time. One source is Ken Clarkson’s code. This is old-style C code which I initially made the mistake of trying to integrate into my own code as a library, eschewing its command line interface. Unfortunately this meant making the code (which is full of global variables) re-entrant; it turns out the code contains subtle crashes and can generate unexpected results if repeatedly run over many inputs. So keeping it as a “black box” and calling it as an exe from my own code seems like the only option – however it does work well in that way, though it is somewhat slower as it has to invoke an exe for every polygon to be created.

Another shortcoming of this particular code is that the alpha shapes produced do not include any “holes”.

For my alpha I’ve chosen 2500, and I’ve also used a multiplier in the command line to reduce the risk of numerical instability with points that are very close to each other (latitude/longitude coordinates going down to 6 or more decimal places).

Another minor issue to note is that when your data includes more than one data item “stacked” at the same spatial point (for example if the data has several points representing the distinct addresses within an apartment building or suchlike), then in the worst case if such stacked points are isolated, you can end up with a cluster that only occupies a single point and hence doesn’t get processed as a polygon.

For very large clusters (bigger than say, 40000 or so) one could potentially also reduce the time for the Alpha Shape calculation by putting a grid over the points and taking only a sample of points in each grid cell.

In the next post we’ll see the results of running the Alpha Shapes code on our clusters. For now, here is the code that calls the Clarkson executable. Note that it contains a Windows-centric function which is not cross-platform, RunExeFromCommandLine.



#include <ogrsf_frmts.h>

#include <fstream>

#include <set>

#include <vector>

#include <iomanip>

#include <windows.h>  // for SHELLEXECUTEINFO etc

#include "node.h"



using namespace std;



bool RunExeWithCommandLine(wstring wstrExelocation, wstring wstrCommand)

{

        SHELLEXECUTEINFO lpExecInfo;

        lpExecInfo.cbSize  = sizeof(SHELLEXECUTEINFO);

        lpExecInfo.lpFile = wstrExelocation.c_str();

        lpExecInfo.fMask = SEE_MASK_DOENVSUBST | SEE_MASK_NOCLOSEPROCESS;     

        lpExecInfo.hwnd = NULL;  

        lpExecInfo.lpVerb = L"open";

        lpExecInfo.lpParameters = wstrCommand.c_str();

        lpExecInfo.lpDirectory = NULL;

        lpExecInfo.nShow = SW_SHOWMINNOACTIVE;  

        lpExecInfo.hInstApp = (HINSTANCE)SE_ERR_DDEFAIL; 

        

        BOOL fSuccess = ShellExecuteEx(&lpExecInfo);

        if (lpExecInfo.hProcess != NULL)

        {

                ::WaitForSingleObject(lpExecInfo.hProcess, INFINITE);

                ::CloseHandle(lpExecInfo.hProcess);

        }

        return fSuccess == TRUE;

}



OGRLinearRing* GetNextRing(vector<pair<int, int> >& rgptsPoly, const vector<Node*>& rgpts, int first, set<int>& new_points)

{

        int rgpt_index = first;

        OGRLinearRing* ring = NULL;



        for (;;)

        {

                int from = -1, to = -1;

                for (vector<pair<int, int> >::iterator it = rgptsPoly.begin(); it != rgptsPoly.end(); it++)

                {

                        if (it->first == rgpt_index)

                        {

                                from = it->first;

                                to = it->second;

                                *it = pair<int, int>(-1,-1);

                                break;

                        }

                        if (it->second == rgpt_index)

                        {

                                from = it->second;

                                to = it->first;

                                *it = pair<int, int>(-1,-1);

                                break;

                        }

                }

                if (from == -1)

                        return NULL;

                OGRPoint point(rgpts[from]->second, rgpts[from]->first);

                if (ring == NULL)

                        ring = new OGRLinearRing;

                ring->addPoint(&point);

                rgpt_index = to;

                new_points.insert(from);

                new_points.insert(to);

                if (to == first)

                        break;

        }

        return ring;

}



int GetFirstNonZero(vector<pair<int, int> >& rgptsPoly)

{

        for (vector<pair<int, int> >::const_iterator it = rgptsPoly.begin(); it != rgptsPoly.end(); it++)

                if (it->first >= 0)

                        return it->first;

        return -1;

}



int GetRings(vector<pair<int, int> >& rgptsPoly, const vector<Node*>& rgpts, vector<OGRLinearRing*>& rings)

{

        set<int> setStartingPoints;

        setStartingPoints.insert(GetFirstNonZero(rgptsPoly));

        for (;;)

        {

                if (GetFirstNonZero(rgptsPoly) < 0)

                        return rings.size();

                set<int> new_points;

                for (set<int>::const_iterator it = setStartingPoints.begin(); it != setStartingPoints.end(); it++)

                {

                        OGRLinearRing* ring = GetNextRing(rgptsPoly, rgpts, *it, new_points);

                        if (ring != NULL && ring->getNumPoints() >= 3)

                                rings.push_back(ring);

                }

                if (new_points.size() == 0)

                {

                        new_points.insert(GetFirstNonZero(rgptsPoly));

                        if (new_points.size() == 0 || *(new_points.begin()) < 0)

                                break;

                }

                setStartingPoints.clear();

                setStartingPoints.insert(new_points.begin(), new_points.end());

        }

        return rings.size();

}



int GetPolygonBoundary(const vector<Node*>& rgpts, int nCluster, vector<OGRPolygon*>& rgNewPolys)

{

        wstring wstrPathForTempFiles = L"C:/temp/";

        wstring wstrTempFileInput = wstrPathForTempFiles + L"sites";

        wstring wstrTempFileOutput = wstrPathForTempFiles + L"sites-alf";

        wstring wstrCommandLineArgs = L"-A -oN -i" + wstrTempFileInput + L" -aa2500 -m10000";



        ofstream os;

        os.open(wstrTempFileInput.c_str(), ios::trunc);

        if (!os.is_open())

                return 0;

        for (vector<Node*>::const_iterator it = rgpts.begin(); it != rgpts.end(); it++)

                os << setprecision(9) << (*it)->second << " " << (*it)->first << endl;

        os.close();



        if (RunExeWithCommandLine(L"C:/Temp/hull/Release/addcl.exe", wstrCommandLineArgs))

        {

                vector<pair<int, int> > rgptsPoly;

                ifstream in;

                in.open(wstrTempFileOutput.c_str());

                if (!in.is_open())

                        return 0;

                char blank[500];

                in.getline(blank, 500); // ignore first line

                while (!in.eof())

                {

                        int n1, n2;

                        in >> n1 >> n2;

                        rgptsPoly.push_back(pair<int, int>(n1, n2));                    

                }

                if (rgptsPoly.size() == 0)

                        return 0;

                vector<OGRLinearRing*> rings;

                if (rgptsPoly.size() <= 2)

                        return 1;

                if (GetRings(rgptsPoly, rgpts, rings) <= 0)

                        return 0;

                int nTotalNumberOfPointsInRings = 0;

                for (vector<OGRLinearRing*>::const_iterator it = rings.begin(); it != rings.end(); it++)

                        nTotalNumberOfPointsInRings += (*it)->getNumPoints();



                // Note if nTotalNumberOfPointsInRings < (int)rgptsPoly.size(), then GetRing only processed 

                //         nTotalNumberOfPointsInRings out of rgptsPoly.size() points.  This is OK, as sometimes

                //         the shape includes a line (polygon with 2 points), which we'll generally ignore.



                for (vector<OGRLinearRing*>::const_iterator it = rings.begin(); it != rings.end(); it++)

                {

                        OGRPolygon* poly = new OGRPolygon;

                        poly->addRing(*it);

                        if (poly->getExteriorRing()->getNumPoints() >= 3)

                                poly->closeRings();

                        rgNewPolys.push_back(poly);

                }

                return rgNewPolys.size();

        }

        return 0;

}

Note that I do not consider the code above to be particularly good; I was still trying to gain an understanding of the output of Clarkson's alpha shape code as I was developing it, and it shows. Anyhow, it has the useful property of working fairly robustly.

The above code is called with code like the following:



vector<OGRPolygon*> rgNewPolys;



map<int, vector<Node*> > mapNodesInEachCluster;



for (vector<Node*>::const_iterator it = vecNodes.begin(); it != vecNodes.end(); it++)

        if ((*it)->m_nCluster >= 0)

                mapNodesInEachCluster[(*it)->m_nCluster].push_back(*it);



for (map<int, vector<Node*> >::const_iterator it = mapNodesInEachCluster.begin(); it != mapNodesInEachCluster.end(); it++)

        if (GetPolygonBoundary(it->second, it->first, rgNewPolys) <= 0)

                cout << "Failure with " << it->second.size() << " points" << endl;



WritePolygonsToMidMifFile("C:/Temp/Polygons", rgNewPolys);

Where the WritePolygonsToMidMifFile function is:



bool WritePolygonsToMidMifFile(string strFileName, vector<OGRPolygon*>& polys)

{

        ofstream osFileMid, osFileMif;



        if (!OpenMidMifAndWritePreamble(strFileName, osFileMid, osFileMif))

                return false;





        int i = 0;

        for (vector<OGRPolygon*>::const_iterator it = polys.begin(); it != polys.end(); it++, i++)



        {

                osFileMid << "i" << endl;

                osFileMif  << "Region  1" << endl

                                << "    " << (*it)->getExteriorRing()->getNumPoints() << endl;          

                for (int n = 0; n < (*it)->getExteriorRing()->getNumPoints(); n++) // first and last point are the same

                {

                        OGRPoint pt;

                        (*it)->getExteriorRing()->getPoint(n, &pt);

                        osFileMif << pt.getX() << " " << pt.getY() << endl;

                }

                int rgb = (40 + ((i * 13) % 215)) * 255 * 255 + (40 + ((i * 4573) % 215)) * 255 + (40 + (i % 215));

                osFileMif  << " Pen (1,2,0) " << endl

                                << "    Brush (" << 2 << "," << rgb << ")" << endl;

        }

        osFileMif.close();

        osFileMid.close();

        return true;

}

We’ll see some of the outputs produced in the next and final post in this series, where we’ll also look at the effect of different values of alpha.

Spatial Clustering in C++ (post 3 of 5) – A Speedier DBSCAN

2011-08-29T12:02:00.002+10:00

Welcome to part 3 of this 5 part blog post on Spatial Clustering using C++. In the two previous posts we read in some map data (Californian census block centroids) and performed some clustering on the points using the DBSCAN algorithm. In this post we’ll investigate how to make the code run faster. This is important if you want this type of algorithm to run in reasonable time over a very large data set (e.g. all of the USA), or you want to run it repeatedly (e.g. to investigate different parameters).

In this investigation I’ll employ some of the techniques for code optimisation that I keep in my mental kit bag, namely: using raw arrays instead of container objects; pre-calculate and cache values wherever possible; use statics to avoid memory overheads; “pool” up frequently re-allocated objects.

But first, the golden rule: always profile. The first thing I’ll do is surround the call to DBScan() with a simple time clocking:



time_t tmStartTime;

time(&tmStartTime);



int nClusters = RunDBScan();



time_t lTimeNow;

time(&lTimeNow);

time_t tsDuration = lTimeNow - tmStartTime;



cout << tsDuration << " seconds to run DBScan" << endl;

The version of DBSCAN presented in the previous post runs in about 86 seconds, so this is the benchmark where we start from (note that is in the Debug build configuration).

Next I do my “poor man’s profiling” – in Visual Studio, I just hit Debug / Break when the main algorithm is running. Breaking a few times gives a picture of where the time is being spent the most. [Incidentally, I love this method of profiling, as it’s just so quick, given that it doesn’t require any external tools or specialised performance builds, yet it still often enables you optimise 90% of the worst bottlenecks in the code].

In this case I find the majority of time is spent allocating memory when pushing items onto the STL vector rgpNeighbourhood2 inside ExpandCluster / GetNodesInRadius.

I tried to use vector’s “reserve” function to allocate memory in larger chunks. However, it’s not clear what the best chunk size to use would be; in fact experimenting with different values can dramatically affect the overall performance. For my test data set, I found the best value was 25, and this reduced the running time to 60 seconds.

Another similar strategy when using STL vectors is to use static and clear; replace the declaration of the vector with the following code:



static vector<Node*> rgpNeighbourhood2;

rgpNeighbourhood2.clear();

vector::clear does not change the vector’s capacity in this version of STL (with MS Visual Studio 8), so the vector is made only once, and it is re-allocated presumably only a handful of times (whenever it's required to get larger than it previously was). This method is better than using reserve, as it doesn’t require any arguments or parameters, and it is faster too – only 53 seconds.

Of course one can also fall back on good old raw C arrays instead of using vectors. I haven’t posted the code for this, but it reduces the time further to 39 seconds. Using raw C arrays does really impact the readability and maintainability of the code and would need to be wrapped up in some nice class wrapper or suchlike in an industrial strength application.

Once this is all done, the bottleneck is now inside kdtree, specifically where it allocates and de-allocates its result objects (objects of type kdres). Looking into the code in this area, I saw the allocation guarded with #defines: USE_LIST_NODE_ALLOCATOR. Looking at the documentation, it advocates switching this tag on. In Windows you will also need to define NO_PTHREADS and comment out the line #ERROR that warns about the code no longer being thread-safe.

Using the list node allocator reduces the time down to an astounding 15 seconds. Finally, running the code in Release instead of Debug reduces the time down to 4 seconds – not bad considering the first naive version took nearly a minute and a half to run! Profiling now shows most time is spent inside kdtree functions to return the nearest nodes, and there’s not much we can do about that.

Summarising what we’ve found:

86 seconds for the first version of the algorithm
60 seconds with reserve of 25 on rgpNeighbourhood2
53 seconds with static instead of reserve
39 seconds with raw C arrays instead of STL vectors
15 seconds with kdtree list_node_allocator switched on (and raw C arrays)
4 seconds in Release configuration (and raw C arrays)

And the lessons learnt (as always, it seems) from optimising code:

always measure timings – bottlenecks are very often not where you’d think they would be
it’s often sufficient to use very simple profiling
it’s frequently surprising how fast you can eventually make the code run

Now that we have our DBSCAN running blazingly fast, we'll in the next post see how to create alpha shapes that will give the output more visual appeal.

Spatial Clustering in C++ (post 2 of 5) – Running DBSCAN

2011-08-26T12:06:00.011+10:00

In the last post in this series I explained how to get and read in some spatial data points for California. In this post we’ll run the DBSCAN algorithm to perform spatial clustering.

DBSCAN is a well-known algorithm that identifies clusters based on density. A major advantage of DBSCAN is that it can identify arbitrary shape objects (ie. it does not presuppose the size, shape, or number of clusters in the data, as some other methods do), and it removes noise during the clustering process. It’s also quite a simple algorithm to implement.

The two parameters controlling DBSCAN are eps and minPts. They work together to define how the clusters should “grow”: if more than minPts points are within a radius of eps of a point, then those neighbouring points are also included in the cluster; each point in the expanded cluster is then subjected to the same test.

Our kdtree data structure provides just what we need for the “neighbours” test with its kd_nearest_range function. However as a side note we should observe that our eps (radius) parameter here is expressed in latitude/longitude distance units – in a real application we would likely want to use a more natural unit such as metres or kilometres. (One can convert, but the conversion is dependent on latitude). One consequence of using lat/long distance, apart from the its non-intuitiveness, is that the radius will be slightly different depending on the latitude of the point.

So, we’ll now show the code to get the points in a given radius – the GetNodesInRadius function:

bool GetNodesInRadius(Node* pt, double dblRadius, int nMinPts, vector<Node*>& rgpNodesFound)

{

        double pos[2] = {pt->GetLongitude(), pt->GetLatitude()};

        kdres* res = kd_nearest_range(kdTree, pos, dblRadius);

        int number_near = kd_res_size(res);

        if (number_near >= nMinPts)

        {

                while (!kd_res_end(res)) 

                {

                        Node* ptNear = (Node*)kd_res_item(res, pos);

                        rgpNodesFound.push_back(ptNear);

                        kd_res_next(res);

                }

        }

        kd_res_free(res);

        return number_near >= nMinPts;

}

And this is called by our DBSCAN code, which is as follows.

void ExpandCluster(vector<Node*>& rgp, int nCluster, double dblEpsilon, int nMinPts)

{

        vector<Node*> rgpNeighbourhood;

        for (int i = 0; i < (int)rgp.size(); i++)

                if (!rgp[i]->m_fVisited)

                {

                        rgp[i]->m_nCluster = nCluster;

                        rgp[i]->m_fVisited = true;

                        rgpNeighbourhood.push_back(rgp[i]);

                }



        for (int i = 0; i < (int)rgpNeighbourhood.size(); i++)

        {

                Node* pNodeNear = rgpNeighbourhood[i];

                vector<Node*> rgpNeighbourhood2;

                if (GetNodesInRadius(pNodeNear, dblEpsilon, nMinPts, rgpNeighbourhood2))

                {

                        // append to rgpNeighbourhood items in rgpNeighbourhood2 that aren't already in rgpNeighbourhood

                        for (int j = 0; j < (int)rgpNeighbourhood2.size(); j++)

                        {

                                Node* pNode = rgpNeighbourhood2[j];

                                if (!pNode->m_fVisited)

                                        pNode->m_fVisited = true;

                                if (pNode->m_nCluster < 0)

                                {

                                        pNode->m_nCluster = nCluster;

                                        rgpNeighbourhood.push_back(pNode);

                                }

                        }

                }

        }

}



int RunDBScan()

{

        double dblEpsilon = 0.0045;

        int nCluster = 0, nMinPts = 8;

        for (vector<Node*>::const_iterator it = vecNodes.begin(); it != vecNodes.end(); it++)

        {

                Node* pNode = *it;

                if (!pNode->m_fVisited)

                {

                        pNode->m_fVisited = true;

                        

                        vector<Node*> rgpNeighbourhood;

                        if (GetNodesInRadius(pNode, dblEpsilon, nMinPts, rgpNeighbourhood))

                        {

                                pNode->m_nCluster = nCluster;

                                pNode->m_fVisited = true;

                                ExpandCluster(rgpNeighbourhood, nCluster, dblEpsilon, nMinPts);

                                nCluster++;

                        }

                }

        }

        return nCluster;

}

Here I’ve implemented a very simple performance improvement over the DBSCAN pseudo-code as it appears in Wikipedia. I don’t add the point P’ to the current neighbourhood N unless it has no cluster already; this guarantees that subsequent iterations over N will only process unvisited points. The results are exactly the same but it runs 2-3 times faster (at least with this STL vector-based implementation).

Even so, this code uses a fairly naive implementation with STL vectors. We’ll see in the next blog post in this series how to refine this implementation and in the process make it run very fast.

With this example data set (California census block data) there tend to be many clusters, including many small ones (though, to be sure, there is at least one truly massive cluster – Los Angeles!). One could obviously cull out the smallest clusters after DBSCAN is run if that is appropriate for your application.

I'm writing out the cluster points to Mid/Mif files and importing them into Tab using MapInfoPro. The code to write the points is:

bool OpenMidMifAndWritePreamble(string strFileName, ofstream& osFileMid, ofstream& osFileMif)

{

        osFileMid.open((strFileName + ".mid").c_str(), ios::trunc);

        if (!osFileMid.is_open())

                return false;

        osFileMif.open((strFileName + ".mif").c_str(), ios::trunc);

        if (!osFileMif.is_open())

                return false;

        osFileMid << setprecision(20);

        osFileMif << "Version 300" << endl;

        osFileMif << "Charset \"Neutral\"" << endl;

        osFileMif << "Delimiter \",\"" << endl;

        osFileMif << "CoordSys Earth Projection 1, 0 Bounds (-1000, -1000) (1000, 1000)" << endl;

        osFileMif << "Columns 1" << endl;

        osFileMif << "  Cluster Integer" << endl;

        osFileMif << "Data" << endl << endl;

        return true;

}



bool WritePoints(string strFileName)

{

        ofstream osFileMid, osFileMif;

        if (!OpenMidMifAndWritePreamble(strFileName, osFileMid, osFileMif))

                return false;



        for (vector<Node*>::const_iterator it = vecNodes.begin(); it != vecNodes.end(); it++)

        {

                osFileMid << (*it)->m_nCluster << endl;

                osFileMif << "Point " << (*it)->GetLongitude() << " " << (*it)->GetLatitude() << endl;

                int c = (*it)->m_nCluster;

                int rgb = (40 + ((c * 13) % 215)) * 255 * 255 + (40 + ((c * 4573) % 215)) * 255 + (40 + (c % 215));

                int size = ((*it)->m_nCluster >= 0 ? 6 : 2);

                osFileMif << "  Symbol (34," << ((*it)->m_nCluster >= 0xFFFFFF ? rgb : size) << "," << 8 << ")" << endl;

        }



        osFileMif.close();

        osFileMid.close();



        return true;

}

Here I've used a cute little one-liner to generate unique colours for each cluster. Non-clustered points ("noise" in DBSCAN terms) are coloured white (with black outline) and drawn smaller than the clustered points.

The results with eps set to 0.0045 and minPts = 8 give 3706 clusters, and 488,761 out of the total of 710,146 points end up in a cluster. Zooming in on the San Francisco area looks like:

Changing eps to 0.009 gives larger clusters, as you would expect. This time there are 3473 clusters, and 592,652 points are clustered. San Francisco looks like:

In the next post I'll show how to significantly speed up DBSCAN, and in further posts demonstrate how to get our output to look like this:

Spatial Clustering in C++ (Part 1 of 5)

2011-08-22T16:29:00.025+10:00

Welcome to the first of a five-part blog post where I work through an end-to-end example of spatial data clustering using C++. I’ll be using a number of tools and techniques, including:

- processing spatial point data using MapInfoPro

- storing points in a kdtree data structure (which will enable fast lookup of nearby points)

- performing spatial clustering using an implementation of the DBSCAN algorithm

- creating Alpha Shapes (think “shrink-wrapped” polygons) to visualise the output

- using OGR for various miscellaneous geometry tasks

Aside from the commercial product MapInfoPro (which is optional for this), all of the tools and codes I use are free/open source with permissive licences. We'll also see (particularly by blog post #3 in this series) that the combination of these tools means our spatial clustering can be blazingly fast to compute.

The simple version of spatial data clustering I'll be doing here is to develop code that will automatically figure out groups (clusters) of densely connected points on a map. A typical application for this might be to automatically find towns or cities from a population map.

The eventual output we'll get looks like the following (and I hope this whets your appetite to follow the entire 5 parts of this blog post):

The above picture shows each cluster found in the data, and coloured uniquely.

The above picture shows the output of running Alpha Shapes to obtain polygons which represent the clusters.

For this exercise I’ll be using a subset of the US Tiger spatial data. In this case I’m choosing California block data. (I did in fact find it very hard to find a freely available file of spatial data points - e.g. addresses for an entire country - to use as my sample data). To get this TIGER data, go to the TIGER download site and Select Block data, California, All Counties.

There are at least two easy ways to get points from this data. The first, and probably quickest (if you have MapInfoPro installed), is to open the block file in MapInfoPro, and go Query / SQL Select, and select CentroidX(obj), CentroidY(obj) from the table, then save this selection as a separate file and export it to MID/MIF format. The MID file will then contain all the centroid points. An alternative shown below is to use OGR-computed centroids. There are 710,146 points in all.

Now let’s write some C++ to read in these points, and store them in a 2-dimensional kdtree, as well as in an STL vector of Node objects. Our kdtree code here is provided by the excellent implementation hosted at Google Code (and consists of just two files, kdtree.h and kdtree.cpp). This code is compiled with Visual Studio but should be relatively cross-platform.

First we define a little helper class to store our points (which we'll call Nodes). The Visited and Cluster members will be used later for DBSCAN.

#include <utility>

class Node : public std::pair<double, double>

{

public:

        Node(double x, double y) : std::pair<double, double>(y,x)

        {

                m_fVisited = false; 

                m_nCluster = -1;

        }

        double GetLongitude() const {return second;}

        double GetLatitude() const {return first;}

        bool m_fVisited;

        int m_nCluster;

};

Next is the main code that reads in the points and stores them in both a kdtree (so that we can get nearest nodes later) and a vector. I've included all the #includes that we're going to need in the code of subsequent blog posts.



#include "kdtree.h"

#include <fstream>

#include <iostream>

#include <iomanip>

#include <vector>

#include <map>

#include <ogrsf_frmts.h>

#include "node.h"



using namespace std;



kdtree* kdTree;

vector<Node*> vecNodes;



bool ReadInPointsFromFile(string strFile)

{

  ifstream in;

  in.open(strFile.c_str());

  if (!in.is_open())

    return false;

  kdTree = kd_create(2);

  while (!in.eof())

  {

    double x, y;

    char comma;

    in >> x >> comma >> y;

    Node* pt = new Node(x, y);

    vecNodes.push_back(pt);

    double xy[2] = {x, y};

    kd_insert(kdTree, xy, pt);

  }

  return true;

}



int main(int argc, char* argv[])

{

  ReadInPointsFromFile("C:/Temp/CentroidPoints.MID");

}

The alternative approach (mentioned above) for getting points from the block data - without needing MapInfoPro - is to read in the shape file data directly and use OGR library functions to calculate the polygon centroids. This code is shown below. You'll need to have OGR installed for this - the easiest way is to install the FWTools binaries. You'll then need to make sure your project includes the FWTools Include file, and links with the appropriate lib.

This approach is a little slower to run than using MapInfoPro. (If you just wanted a fast approximation of the block centroids, however, you could probably just take the centroid of the bounding box of the polygons as they are read in.)



bool ReadInPolygonsFromFile(string strFile)

{

  kdTree = kd_create(2);

  OGRRegisterAll();

  OGRDataSource* poDS = OGRSFDriverRegistrar::Open(strFile.c_str(), false);

  if (poDS == NULL || poDS->GetLayerCount() <= 0)

    return false;

  OGRLayer* poLayer = poDS->GetLayer(0);

  OGRFeature* poFeature;

  while ((poFeature = poLayer->GetNextFeature()) != NULL)

  {

    OGRGeometry* poGeometry = poFeature->GetGeometryRef();

    if (poGeometry != NULL)

    {

      OGRwkbGeometryType nType = wkbFlatten(poGeometry->getGeometryType());

      OGRFeatureDefn* poFDefn = poLayer->GetLayerDefn();

      long lType = wkbFlatten(nType);

      if (lType == wkbPolygon)

      {

        if (static_cast<OGRPolygon*>(poGeometry)->get_Area() != 0)

        {

          OGRPoint pt;

          static_cast<OGRPolygon*>(poGeometry)->Centroid(&pt);

          Node* node = new Node(pt.getX(), pt.getY());

          vecNodes.push_back(node);

          double xy[2] = {pt.getX(), pt.getY()};

          kd_insert(kdTree, xy, node);

        }

      }

      else if (lType == wkbMultiPolygon)

      {

        OGRMultiPolygon* pmpg = static_cast<OGRMultiPolygon*>(poGeometry->clone());

        for (int nPoly = 0; nPoly < pmpg->getNumGeometries(); nPoly++)

        {

          OGRGeometry* poSubGeometry = pmpg->getGeometryRef(nPoly);

          if (poSubGeometry != NULL && static_cast<OGRPolygon*>(poSubGeometry) != NULL)

          {

            if (static_cast<OGRPolygon*>(poSubGeometry)->get_Area() != 0)

            {

              OGRPoint pt;

              static_cast<OGRPolygon*>(poSubGeometry)->Centroid(&pt);

              Node* node = new Node(pt.getX(), pt.getY());

              vecNodes.push_back(node);

              double xy[2] = {pt.getX(), pt.getY()};

              kd_insert(kdTree, xy, node);

            }

          }

        }

      }

      else  //  maybe lType == wkbPoint or wkbLineString ?

        return false;

    }  

    OGRFeature::DestroyFeature(poFeature);

  }

  return false;

}

int main(int argc, char* argv[])

{

  ReadInPolygonsFromFile("C:/Temp/tl_2010_06_tabblock10/tl_2010_06_tabblock10.shp");

}

Stay tuned for the next post, where we’ll run the DBSCAN algorithm and start to see some nice outputs!

js/css resource serving in python apps with Fanstatic

2011-07-27T16:33:00.005+10:00

I've just been checking out Fanstatic, a resource publishing/static file serving solution for wsgi python apps. I've been contemplating something like this as our javascript and css dependencies are getting more complex. It would also be useful to have some form of automatic cache invalidation so users don't have to do a special browser refresh when we update our applications.

It's easy to set up with CherryPy

from fanstatic import Fanstatic

if __name__ == "__main__":

    app = cherrypy.Application(Root())
    app.wsgiapp.pipeline.append(('repoze.who', setup_auth))
    app.wsgiapp.pipeline.append(('beaker', setup_session_storage))
    app.wsgiapp.pipeline.append(('fanstatic', Fanstatic))
    cherrypy.quickstart(app, config='workbench.conf')

after that is done, you can jquery.need() in the widget/template that needs jquery, and similar for our other dependencies. Has anyone else used fantastic? What are other solutions to dependencies and serving js/css? Is there an easier and better solution? Wrapping new libraries for fanstatic looks like a bit of effort but I haven't explored it much yet.

Loki

Multi-server SSL with CherryPy

2011-07-08T15:22:00.002+10:00

It's quite hard to see in the CherryPy docs how to support http and https from a single app. It turns out combining a few of the posts on the tools wiki is the way to go.

Starting the app looks like:

if __name__ == "__main__":

    from cherrypy import _cpserver
    from cherrypy import _cpwsgi_server

    secure_server = _cpwsgi_server.CPWSGIServer()
    cherrypy.config.namespaces['secure_server'] = lambda k, v: setattr(secure_server, k, v)
    secure_server.bind_addr = ('0.0.0.0', cherrypy.config['secure_server.socket_port'])

    adapter = _cpserver.ServerAdapter(cherrypy.engine, secure_server, secure_server.bind_addr)
    adapter.subscribe()

    app = cherrypy.Application(Root())
    app.wsgiapp.pipeline.append(('repoze.who', setup_auth))
    app.wsgiapp.pipeline.append(('beaker', setup_session_storage))
    cherrypy.quickstart(app, config='workbench.conf')

And the important config file options are:

[global]
server.socket_host = "0.0.0.0"
server.socket_port = 8080
server.thread_pool = 10
secure_server.socket_port = 8443
secure_server.ssl_certificate = "certs/self_sign_certificate.crt"
secure_server.ssl_private_key = "certs/self_sign.key"

Loki

Biarri provides Network Design Optimisation to NBN Co.

2011-06-23T14:49:00.000+10:00

NBN Co today announced that they awarded an 8 year contract to Biarri to supply network design optimisation to support development of efficient, lower cost network construction plans.

Biarri is to provide optimisation technology to NBN Co for fibre network design optimisation. The optimisation engine quickly generates low-cost fibre network designs based on the requirements of the reference architecture. It can determine optimal fibre area boundaries, the position of fibre hubs, and the layout and route of distribution and local fibre.

See the NBN Co site for more detail.

Logic of the Fibre Optic Network Design Tool

Optimisation: Striking the Right Balance

2011-06-03T09:59:00.006+10:00

One of the guiding principles we use in commercial mathematics is to "Model Conservatively, but Optimise Aggressively". This means that the problem domain should be modeled with sufficient "fat" in the data to ensure that the results are both legal and robust; but given this, we should then seek to apply the best (fastest and highest quality) solution approach that we can get our hands on.

Optimising aggressively can sometimes have it's downfalls, though, if taken too literally. I've been doing a few experiments with numerical weightings of the objective function in a Vehicle Routing problem, where this issue is readily apparent. (Actually it is a Vehicle Routing problem with time windows, heterogeneous fleet, travel times with peak hours, both volume and weight capacities, and various other side constraints).

Our Vehicle Routing uses travel times (based on shortest paths through the street network) that are characterised by distance and duration. Durations can vary due to different road speeds on different types of streets (highways vs suburban roads for example). This leads to the question of how (on what basis) to optimise the vehicle routes - given that the optimisation has already to some extent minimised the number of vehicles and created well-clustered routes - what is the most desirable outcome for KPIs in terms of duration and distance?

In one experiment I've tried three different weightings for the duration (cost per hour) while keeping the cost per distance constant. I've run three values for this cost per hour - low, medium, and high weightings - on real-life delivery problems across two different Australian metropolitan regions.

		Region 1
	Total Duration	Driving Duration	Distance
Cost/hour
Low	74:47	24:38	708
Medium	72:45	23:55	712
High	72:58	23:42	768

		Region 2
	Total Duration	Driving Duration	Distance
Cost/hour
Low	113:54	46:44	1465
Medium	107:51	41:36	1479
High	108:51	43:49	1518

From these results, there is a (more-or-less) general correspondence between distance and the driver cost per hour as you would expect. However, if you push one weighting too far (ie. optimise too aggressively or naively), it will sometimes be to the detriment of all the KPIs as the optimisation will be pushing too strongly in one direction (perhaps it is outside the parameter space for which it was originally tuned, or perhaps it pushes the metaheuristic into search-space regions which are more difficult to escape from). This is most acutely seen in Region 2 when using the high cost per hour value. Conversely if you drop the cost per hour to a low value, the (very modest) reduction you get in distance is very badly paid for in terms of much longer durations. What is most likely happening in this case is that the routes are including much more waiting time (waiting at delivery points for the time windows to "open"), in order to avoid even a short trip (incurring distance) to a nearby delivery point that could be done instead of waiting.

The problem of striking the right balance is most acute with metaheuristics which can only really be tuned and investigated by being run many times across multiple data sets, in order to get a feel for how the solution "cost curve" looks in response to different input weightings. In our example, an in-between value for cost per hour seems to strike the best balance to produce the overall most desirable KPI outcome.

Successfully Implementing Commercial Mathematics

2011-05-27T15:52:00.000+10:00

We have grown Biarri based on successfully bringing the power of mathematics to bear on business. We try hard to make operations research easier for business to digest. Over the last few years there are a few things we have learned that are worth pointing out.

1. Target businesses that need analytical support and know they need it. Many businesses now have lots of data but that they lack the analytical ability to utilise that data. The businesses we target understand that smart decisions are based on 'results gleaned from algorithmic insight and executed with the confidence that comes from really doing the math'. 'Analytics' is now part of the business improvement conversation however you need to target businesses that have complex problems that are solvable and need solving. Additionally, we have always had better success when the problem has a clear owner who will benefit directly from solving the problem.

2. Differentiate your offering. 'Optimisation' and now 'Analytics' are overused terms. As applied mathematicians we can certainly bring some real science to the table. But this is not enough. What has emerged as very important is our ability to make the mathematics digestible. We have tried hard to deliver 'Accessible Optimisation'. We believe optimisation is only powerful if it is:

Easy to apply to real world operations
Easy to understand (intuitive)
Easy to access (usable)
Affordable with minimum capital spend
Fast and reliable

3. Focus on implementation practicalities. You need to keep your eye on the basics of project management. Some of the key issues for us have been:

Work hard to keep the scope tight by focusing on the 20% of features that deliver 80% of the value.
Use simple frameworks that everyone can understand. All our custom models and tools use a simple linear workflow/logic around a smart engine.
Stay close to your customer and practice the art of no surprises.
Excel is useful for modelling. It is available on every desk top, flexible, familiar, easy to use and accepted by non-technical managers.
Develop and implement prototype models quickly and have a simple way of upgrading them when needed. For example, we always keep our prototype engines separate to the front-end so we can easily port them to different solution frameworks.

Custom Excel Charts: Scatter Plots

2011-05-16T10:32:00.004+10:00

I am often required to create some very specific charts for clients that are far outside of the capability of regular Excel charting functionality. I use VBA to create these custom charts in a dynamic way to visualise outputs of mathematical modelling. These outputs are important to clients as they allow them to view the complex outputs in simple, familiar ways.

I have in the past created L-Shaped Graphs for clients that helped optimise the use of space for a dynamic Excel Dashboard that needed to be constrained to an A4 print out.

More recently, I have been required to create a custom ‘March Chart’ that plots the movements of different work crews over space and time. I received a print out of an old ‘March Chart’ that had been manually created by the client, and was asked to replicate it exactly in Excel so that it could be viewed as an output of a mathematical model.

THE POWER OF SCATTER PLOTS

When creating a custom chart in Excel, you must use one of the existing charts as a template. This template can then be twisted and abused into doing what you need through tricks and manipulation within VBA.

For the March Chart project, I discovered just how very useful a scatter plot can be as the template for an advanced Excel chart. The scatter plot can be thought of as a blank canvas on to which x,y data points can be plotted. Each data point can then be joined by a line to another data point, creating a straight line. A series of sequential data points with straight lines joining them called a ‘data series’ can be used to create paths and even ‘curves’. With this functionality in mind, you essentially now have a blank canvas that you can draw anything on in Excel with VBA!

This all seems far too good, so Excel had to put some barriers in our way to make it a little trickier.

For starters, there is the ever present issue of backward compatibility to Excel 2003. The object model used to create charts was changed from Excel 2003 to Excel 2007, so it is very important to be aware of this, and find out early if your client is using Excel 2003. Fortunately for myself, the client for the ‘March Chart’ had upgraded to at least Excel 2007 so this was not an issue for me.

More specific to the advanced use of scatter plots, Excel limits you to 255 data series on a chart. Each sequence of coordinates that join together to make a continuous line are defined as a data series. This means that you can only have 255 distinct continuous lines on your chart.

Each of these data series, can only consist of 255 data points. Unless your chart is extremely big and going to be shown on a very large screen, this is pretty easy to manage. There are two options:

o Split the data series into multiple data series. If you make the two data series the same colour, and begin the second data series at the last point of the first data series then they will appear to be a continuous line. This is a good option if you are not close to the limit of data series, however can create more issues to overcome as an additional legend entry will be created.

o For any data series that has more than 255 potential data points, skip over every second, third, etc. depending on the actual number so that the number of data points implemented to the actual plotted data series remains under 255. This is my preferred method as no additional work is required in fixing the legend and is in my opinion neater.

The following is an example output of the March Chart showing different colours for different work crew types. You can track the geographic location over time of each work crew very simply with this visualisation.

Excel 2D Data Lookup using INDEX with MATCH

2011-05-16T09:42:00.006+10:00

INDEX with MATCH enables lookup of a cell in with an unknown reference in a matrix with row and column labels.

The two formulas:

INDEX

Syntax: INDEX(array, row_num, [column_num])

Index is very simple and the syntax is fairly self explanatory.

INDEX takes an array of data as the first parameter, a row number as the second and a column number as the third. The function then return the value in the array referenced by the row and column specified.

An example array in Excel:

To return ‘Black’ using the INDEX function use the formula:

=INDEX(B2:D4,3,2)

Note that the second parameter refers to the 3^rd row in the array B2:D4 and not the 3^rd row of the spreadsheet. This applies for the column reference also.

When used like this, INDEX appears to be fairly limited in its practical applications, as if you already know the row and column number of a cell in the array, you can generally reference the cell directly.

However, when used with MATCH, INDEX becomes one of the most powerful and useful functions in Excel.

MATCH

Syntax: MATCH(lookup_value, lookup_array, [match_type])

MATCH is again quite simple on its own however requires a little more consideration and care than INDEX.

MATCH takes a lookup value as its first parameter and a lookup array as the second. For our purposes the lookup array needs to be a single dimension. This can be row wise or column wise.

MATCH then returns the index of the lookup value if it is found in the lookup array.

However false positives are easy, as the third parameter of MATCH is the [match_type] which takes one of three values:

The default value is 1, and this can cause issues as if you forget to add this optional parameter then you could very easily get incorrect lookups occurring that are hard to notice. This is because the lookup will still return a value even if it can not find what you are looking for. MATCH will return some value that is ‘Less than’ the lookup value, which is very rarely what you actually want. It is important when using MATCH with INDEX for matrix lookups to remember to set match_type parameter equal to 0.

An example array in Excel:

To return the row index of the array for the value “Colour 2”:

=MATCH("Colour 2",A2:A4,0)

This will return: 2

Note that this is the row number of the lookup value “Colour 2” in the lookup array and not in the worksheet.

INDEX with MATCH

MATCH is used to provide logic to the row_num and column_num parameters of INDEX.

If a matrix in Excel has labels for both the rows and columns, then INDEX with MATCH can lookup a specific cell in the matrix with ease.

An example of such a matrix:

INDEX with MATCH can be used with this matrix to answer questions such as ‘What is Colour 2 of Group 3?

Repeating from above, the syntax for INDEX is:

INDEX(array, row_num, [column_num])

The input array for this example is B2:D4.

The row_num and column_num parameters need to be found using MATCH as shown above.

row_num: MATCH("Colour 2",A2:A4,0) will return 2

col_num: MATCH("Group 3",B1:D1,0) will return 3

These can then be put into INDEX to create our formula:

=INDEX(B2:D4,MATCH("Colour 2",A2:A4,0),MATCH("Group 3",B1:D1,0))

This will return pink.

I now use INDEX with MATCH nearly every time I do data analysis and manipulation in Excel. Clients often have data in 2D Matrices, and it is no longer a hurdle for me to use this data with ease.

Editing vector layers with Quantum GIS

2011-05-05T10:27:00.004+10:00

I've been looking around for an open source alternative to the excellent Desktop GIS Mapping Application MapInfo Pro. I've installed and played with both Quantum GIS (QGIS) and MapWindow. The latter is a little bare bones and does not seem to include any geometry editing, so I've been focussing instead on QGIS.

QGIS is a great application, though it is quite noticably slower to render very large mapping layers than MapInfo Pro. However - and this is, strangely, not very well known or documented - QGIS will only let you edit layers which are Shapefiles. There is no documentation that I can find which says why the "Toggle Editing" function is always disabled for other files - this is very confusing/frustrating if you don't know about this Shapefile limitation.

There is a converter within QGIS (which is really just a UI on ogr2ogr) to convert between TAB and SHP formats, but Shapefile layers are rather limited because you can't mix geometries within a layer - both nodes (points), and lines, for example. This is rather a problem because the map data I have has both points and lines - the points are needed to style the line endpoints and allow for data attributes on them. Separating points and lines into different layers in this context is a bit of a pain. Shapefiles also have various other limitations - for example, field names can be no longer than 10 characters.

Perhaps hooking up QGIS to an SQLite database with Spatialite extension might make managing the map data layers more streamlined - one for a rainy day...

Installing Cplex in Linux Ubuntu

2011-05-05T09:38:00.002+10:00

I've been trying to install ILOG's Cplex product on a 64 bit Ubuntu machine. There are a few small hiccups I encountered.

First, try to install as root with

./cplex_studio122.linux-x86.bin

On my machine I'm installing to /opt/ILOG/CPLEX_Studio122. If you get an error like "jre/bin/java: not found" then you need the "32 bit libs" package:

apt-get install ia32-libs

(You may also need to set the path with LD_LIBRARY_PATH=/usr/lib32). The 32 bit libraries seem to be required only for the installer (at least, they are not needed for programs that just link with the static Cplex libs).

After installing, you may get a build error (running as a non-root user) that it can't find the cplex header files. Try ls /opt/ILOG/CPLEX_Studio122 and see if there are permission denied messages - the installation seems to screw up some permissions on this folder, however this is easily fixed with chmod +r /opt/ILOG/CPLEX_Studio122.

If you are using COIN-OR's Osi class OsiCpxSolverInterface you will also need the following at the top of your OsiCpxSolverInterface.cpp file:

#include "/opt/ILOG/CPLEX_Studio122/cplex/include/ilcplex/cplex.h"

A typical Makefile snippet which includes Cplex and COIN/Osi might then look like:

my_objects = YourFile.o OsiCbcSolverInterface.o OsiCpxSolverInterface.o

CPPFLAGS = -fPIC -I/usr/include -I/usr/include/coin -DNDEBUG -I/opt/ILOG/CPLEX_Studio122/cplex/include/ilcplex

LIBFLAGS = -l:libCbc.so -l:libCbcSolver.so -l:libCoinUtils.so -l:libOsi.so -l:libOsiClp.so -l:libClp.so -l:/opt/ILOG/CPLEX_Studio122/cplex/lib/x86-64_sles10_4.1/static_pic/libcplex.a -l:/opt/ILOG/CPLEX_Studio122/cplex/lib/x86-64_sles10_4.1/static_pic/libilocplex.a

my_program : $(my_objects)
g++ -Wall -fPIC -shared -o my_program $(my_objects)

.PHONY : clean

clean :
rm *.o my_program

Function within a Function in C++ - Everything In Its Right Place

2011-03-18T11:49:00.011+11:00

I have been coding in C++ for a number of years. Only recently did I find out that there is an indirect way to write a function within a function. I have found myself using this little workaround more and more recently. When coupled with the use of a static variable within a function, the ability to have "Everything In Its Right Place" i.e. avoid unnecessary functions and variables in the global namespace, is particularly appealing to my inner code nazi.

The below code snippet shows an example...

Note

- I have only tested the code snippet below for compilation using Visual Studio

- Visual Studio's debugger does not allow you to "watch" spndI when execution is inside the function NodeIsOfSameTypeAsIButNotI

- I have formatted the code for the narrow width afforded this blog! So the example looks more verbose than it would otherwise be, and the blogger may mangle some of the code when it is posted (causing me to re edit the post over and over until google gets it right)

struct CNode

{

// shell struct to enable code compilation

double GetLongitude() const {return 0;}

double GetLatitude() const {return 0;}

double GetType() const {return 0;}

};

template class CSomeClassHoldingNodesForFastFinding

{

public:

typedef bool(*NodeCanBeClosest)(const T*);

template inline T* FindClosest

(

double dblLongitude,

double dblLatitude,

NodeCanBeClosest pNodeCanBeClosest

)

{

return NULL; // Really return a T* not Null

}

};

void FindClosestNodePairsAndDoSomethingWithThem

(

vector& rgpnd,

CSomeClassHoldingNodesForFastFinding* pFastFindNodes

)

{

static CNode* spndI = NULL;

for (size_t nNodeI = 0; nNodeI != rgpnd.size(); nNodeI++)

{

spndI = rgpnd[nNodeI];

struct Func

{

static bool NodeIsOfSameTypeAsIButNotI(const CNode* pnd)

{

return

pnd != spndI

pnd->GetType() == spndI->GetType();

}

};

CNode* pndJ = pFastFindNodes->FindClosest

(

spndI->GetLongitude(),

spndI->GetLatitude(),

&Func::NodeIsOfSameTypeAsIButNotI

);

// .. Do something with spndI and pndJ

}

The business of Biarri

2011-02-08T21:56:00.002+11:00

At Biarri we build mathematical optimisation products and models and we have also been building up the business over the past two years from very humble beginnings. Below are a few more of the business lessons we have learnt growing a tech start up in Australia.

1. Get some rhythm – we have used the Rockefeller Habits checklist from Verne Harnish – it gives each day, week, month, quarter and year a target and keeps the team accountable for key actions. It also means we are always watching the cash burn. To quote Verne “cash is oxygen”. This is important as we were also told very early that costs have a way of taking care of themselves if not very closely watched – this is true.

2. Meet your competition – how else will you really know how to differentiate and where their weaknesses are?

3. Give people space to deliver – give good people the space to deliver and they almost always will. People don’t need rules or process – they need a challenges and support. We have been amazed at the great work our super smart university students produce on our projects. Our customers also love it.

4. Keep low cost – we have had many big corporate guys come to our offices and comment on the poor paint job, old carpet or our back alley address. They don’t seem to get it. You don’t need to be on level 24 to produce a great product and you can certainly sell it cheaper if you are not.

5. Stay close to customers – make the calls every day to stay close to people and opportunities arise. We find that you need to be out speaking to people, catching up for coffee and on the phone to grow a business like ours.

Using libCurl to make web requests from a C++ program

2010-12-30T15:46:00.005+11:00

This blog post is to share how to make a C++ program using the libcurl C interface under Windows/Visual Studio, as it is a little non-trivial.

The libcurl API lets you transfer files to a server and call exposed functions outside of a browser context. Because we are using CherryPy with basic user authentication as our web platform, I needed to be able to call various functions such as the one I'm using in this example: "get_result". The equivalent Linux curl command line would be:


curl --user user:password  "httpss://www.insertyourURLhere.com/get_result?id=30417&task=9332"

where obviously you'd change "insertyourURLhere", use the appropriate user name and password, and I've assumed "get_result" takes two arguments, "id" and "task".

On Windows with Visual Studio, after downloading the libcurl source, you need to first build the libcurl lib. To achieve this, you need to change the C++ Preprocessor definitions by removing "_USRDLL", and adding instead "CURL_STATICLIB".

Now, in your own program, add the C++ preprocessor definition "CURL_STATICLIB", link with libcurl.lib as well as with ws2_32.lib, winmm.lib and wldap32.lib. This will avoid link errors.

The following C++ code shows how to make requests to the web server "get_result" function, using HTTP POST. It demonstrates using a callback function (to retrieve the text of the server's response), passing the URL encoded arguments (encoded using curl_easy_escape), authenticating using basic authentication, and gathering timing diagnostics.


#include "curl.h"

using namespace std;

#define YOUR_URL "https://www.insertyourURLhere.com/"
#define USER_AND_PWD "user:password"

static string gs_strLastResponse;

// Callback to gather the response from the server.  Comes in chunks (typically 16384 characters at a time), so needs to be stitched together.
size_t function_pt(void *ptr, size_t size, size_t nmemb, void * /*stream*/)
{
 gs_strLastResponse += (const char*)ptr;
 return size * nmemb;
}

bool CallServerWithCurl(string strData1, strData2, string& strErrorDescription)
{
 CURL* curl = curl_easy_init();
 if (curl == NULL)
 {
     strErrorDescription = "Unable to initialise Curl";
     return false;
 }

 curl_easy_setopt(curl, CURLOPT_URL, (YOUR_URL + "get_result").c_str());
 curl_easy_setopt(curl, CURLOPT_HTTPAUTH, CURLAUTH_BASIC);
 curl_easy_setopt(curl, CURLOPT_USERPWD, MY_USER_AND_PWD);    // set user name and password for the authentication

 char* data1 = curl_easy_escape(curl, strData1.c_str(), 0);
 char* data2 = curl_easy_escape(curl, strData2.c_str(), 0);
 string strArguments = "id=" + data1 + "&task=" + data2;
 const char* my_data = strArguments.c_str();

 curl_easy_setopt(curl, CURLOPT_POSTFIELDS, my_data);
 curl_easy_setopt(curl, CURLOPT_POSTFIELDSIZE, (long)strlen(my_data));   // if we don't provide POSTFIELDSIZE, libcurl will strlen() by itself

 curl_easy_setopt(curl, CURLOPT_VERBOSE, 1L);     // enable verbose for easier tracing

  gs_strLastResponse = "";
  curl_easy_setopt(curl, CURLOPT_WRITEFUNCTION, function_pt);        // set a callback to capture the server's response

 CURLcode res = curl_easy_perform(curl);

 // we have to call twice, first call authenticates, second call does the work
 res = curl_easy_perform(curl);

 if (res != CURLE_OK)
 {
     strErrorDescription = "Curl call to server failed";
     return false;
 }
 if (!DoSomethingWithServerResponse(gs_strLastResponse))
 {
     strErrorDescription = "Curl call to server returned an unexpected response";
     return false;
 }

 // extract some transfer info
 double speed_upload, total_time;
 curl_easy_getinfo(curl, CURLINFO_SPEED_UPLOAD, &speed_upload);
 curl_easy_getinfo(curl, CURLINFO_TOTAL_TIME, &total_time);
 fprintf(stderr, "Speed: %.3f bytes/sec during %.3f seconds\n", speed_upload, total_time);

 curl_easy_cleanup(curl);

 return true;
}

Libcurl building and linking tips were found at http://curl.haxx.se/mail/lib-2007-04/0120.html and elsewhere.

Computer Science Journal Recommendations

2010-12-23T16:17:00.005+11:00

I was just contemplating this last night and I think my slightly before new years resolution is to read one journal paper a week as part of my "Biarri Day learning." :)

I'm going to start with The Journal of AI Research and Kybernetika because they publish everything for free online and are reputable / well ranked. Well... and maybe because Kybernetika has a fantastic sounding name. :)

If anyone else can recommend a useful and readable Computer Science / AI related journal please comment!

Loki

MapInfo tips

2010-12-23T12:02:00.002+11:00

I've had reason to manipulate some spatial data recently using MapInfo Professional, and picked up a few tips that I thought might be useful to someone else one day:

Some good tips for using SQL Select can be found here: http://www.sgsi.com/MIUserGroup/SGSI_SQLSelect.htm.

You can use SQL conditions that select objects by their spatial properties. Use syntax like: Str$(obj) = "line", and Str$(obj) = "point", Str$(obj) = "region" etc.

To change the projection, use File / Save Copy As and click the Projection button.

Get used to using the arrow keys to pan around the map, and the mousewheel to zoom in and out - this saves having to switch to the Grabber tool all the time to move around.

Press "S" to go into "Snap" mode. This is very useful when creating lines, to make sure they attach to nodes (note the nodes don't have to be in the same layer when your lines snap). When in Snap mode you will see the text "SNAP" in the status bar at the bottom of the MapInfo main window.

The menu item Map / Change View is useful to go straight to a given lat/long coordinate. (It centres the view on your entered coordinate). You will need to be careful however if you have multiple layers open that have different coordinates or projections.

To merge two different layers:

* First, the tables have to have the same structure. Use Table / Maintenance / Table Structure to add and remove fields as required.
* Then use either Table / Update Column or Table / Append Rows To Table. If using Update Column on node/arc layers (to copy data from the node or arc objects of one table to the corresponding objects on another), I found using Join with Intersects worked best.

To separate a single layer that contains nodes and arcs into two layers:

* Run an SQL Select query with the condition Str$(obj) = "line" - this will select only the arcs.
* Then Save the selection to a separate file (which will be the arcs layer), and delete the selected arcs from your original layer (leaving only the nodes).

Database vs pure python performance

2010-12-07T12:14:00.003+11:00

I was just looking through some of our common code and noticed some maths code in python that I knew could be replaced by a function in our geo-database, Postgis. However I wondered what the performance difference between calling out to some C code via the database and doing a simple bit of maths in python would be. The code has no loops, is simple and should be very fast.

It turns out the database approach is around 100X slower...



>>> timeit.Timer('common.lat_long_distance_db(db_wrap, 144.865906, 144.865237, -37.836183, -37.831359 random.random())', 
         "from __main__ import common, db_wrap, random").timeit(10000)
2.001802921295166
>>> timeit.Timer('common.lat_long_distance(144.865906, 144.865237, -37.836183, -37.831359 random.random())', 
         "from __main__ import common, random, db_wrap").timeit(10000)
0.02936482429504394

I would suggest this isn't the right way to use the DB! :)

If we write a function where we just select a constant value from the database we can find the call over head of the database calls and remove the postgis process and parsing cost.



timeit.Timer('lat_long_distance_db_just_select(db_wrap, 144.865906, 144.865237, -37.836183, -37.831359 random.random())', 
         "from __main__ import db_wrap, random, lat_long_distance_db_just_select").timeit(10000)
1.4460480213165283

Large percentage of overhead.

Seems to suggest that I should look into the difference between selecting out the data, apply the python function then updating the db with results. I suspect apply the postgis functions to data that is already in the db instead of looping in python will be best, but empirical data is enlightening!

Cherrypy and coping with Internet Exporer Cache issues

2010-12-03T16:39:00.002+11:00

Ahh, IE caches Ajax requests. This is rather confusing when you call a function and expect it to execute something on the server but instead it immediately returns. Luckily, it's very simple to fix when using cherrypy. Just add the following code to your .conf file for the functions you don't want to cache, or apply it globally.



tools.response_headers.on = True
tools.response_headers.headers = [('Expires', 'Sun, 19 Nov 1985 05:00:00 GMT'),
('Cache-Control', 'no-store, no-cache, must-revalidate, post-check=0, pre-check=0'),
('Pragma', 'no-cache')]

Loki

Array Comparison in PostgreSQL

2010-12-03T16:12:00.003+11:00

I just had an input table, where I wanted to map multiple items to a single code.

I had foo and bar which I wanted to map to baz. They could be in any order though, (bar and foo) so I couldn't just aggregate a string and compare those.

Instead I used array_agg and array_sort



select i.id, m.type from 
    (select id_mot1 AS id, array_sort(array_agg(annoying_types)) AS input_type, 
        from 
        (select m.* from import_temp as m order by m.id_input, m.id asc) AS t 
        group by id_input) as i, 
type_map as m where i.input_type = array_sort(m.input_type) ;

I create an array of the grouped values using array_agg, sort it using array_sort and then compare it to the sorted mapping. Sorting on insert into the mapping would make more sense though.

array_sort is from the post by David Fetter at PostgreSQL SQL Tricks.

Loki

Some Linux IDEs for C++

2010-11-18T18:10:00.006+11:00

I recently went on a short excursion looking for good, lightweight IDEs for C++ in Linux, similar to my adventures with Python IDEs. My (modest) needs are similar to my Python IDE needs, specifically however I want very good step-through debugging support, and in terms of a lightweight solution, would like a tool that doesn't have a big project/workspace overhead with the IDE creating many of it's own files (our C++ programs all have simple makefiles, with debugging info turned on and code optimisation suppressed by compiling with the -g and -O0 flags).

Here is what I found:

GDB is the GNU debugger used by GCC and actually has a text interface for step-through debugging that is actually usable - if you are very patient. There is also a graphical interface called DDD, however I did not try this.
Anjuta DevStudio is an IDE for GNOME. Debugging is activated with a plugin found under Edit/Preferences/Installed Plugins. This IDE worked OK for me in 32 bit Linux, but failed (shutting itself down) when trying to debug under 64 bit Linux.
CodeLite is a cross-platform IDE, that seemed to work well at first, however I could not get it to respect my program's command line arguments.
Geany is an IDE with a clean presentation that I have used in the past just to look at C++ code; debugging is supposedly achieved with a plug-in called GeanyGDB, however I did not try it.
Code::Blocks is a very mature IDE that is a little chunky-looking in appearance (though it's download size suggests it is more lightweight than the IDEs listed above), and seemingly requires you to make various project files (at least that seemed the easiest way to get it to a point where I could use it for debugging). However, the debugging does work rather well (and it even supports the Mighty Middle Click). Though it does not seem to show variable internals for more complex variables (dereferencing STL iterators would be nice, but maybe I'm missing something there), this will be my IDE of choice in Linux for a little while at least, I think.

Trying out these IDEs in Linux has, if nothing else, emphasised to me just how impressive, powerful and usable Microsoft's Visual Studio is (for Windows only, of course). I have used Visual Studio in it's various incarnations since 1996, and it keeps upping the bar even with the 2010 version with awesome debugging features like "Pin to source", which lets you make little variable watch windows right alongside the line of code that it is relevant for. It offers so much at a glance (local variable values, current stack, etc) in an elegant way, and one quickly becomes used to being able to look deep into the internals of STL data types, change the line of execution or a variable's value while debugging, and the like.

One thing that would be really good is if all these IDEs standardised their debugging hotkeys; as the table below shows (where I have also included Eric and Firebug), they are generally all different. I think Eric's scheme is the best, as it seems logical and doesn't require any contortions with Shift and Control keys.

	Start	Continue	Step	Step into	Step out	Stop
Anjuta	F4	F4	F6	F5	Shift-F5	none
Codeblocks	F8	Ctrl-F7	F7	Shift-F7	Ctrl-Shift-F7	none
Visual Studio and CodeLite	F5	F5	F10	F11	Shift-F11	Shift-F5
Eric	F5	F6	F7	F8	F9	F10
Firebug	n/a	F8	F10	F11	Shift-F11	n/a

Tech Start Up Business Lessons

2010-10-11T16:22:00.004+11:00

Biarri commenced as a commercial maths start up almost two years ago. In this time we have learnt a lot.
We focus hard on deploying optimisation and quantitative analytics in accessible ways to deliver the power of mathematics quickly and cheaply. Our just launched workbench (www.biarriworkbench.com) is a great example of this – providing monthly low cost rental of powerful maths engines available over a browser.

While we have been building products and models we have also been building a business and have learnt a few things along the way. Below are a few of the business lessons we have learnt growing a tech start up in Australia.

Be in the cloud – because we were delivering our optimisation workbench using the cloud, we sought out cloud services for our internal business needs. Our accounting software, CRM, email and timesheets are all rented from Software as a Service companies. We learnt a lot about what makes a good web app by using these services and we saved a lot of capital cost upfront. Specifically let me say that SAASU (www.saasu.com.au) is a really good accounting system for a small business – much easier to use than MYOB or quicken in my view.
Always push back against one-sided contract terms from big corporates – we find almost always you will get at least what you ask for. In house lawyers and legal departments will always try it on, especially when they are dealing with a small business – push back hard there is always some flex
Not all phone companies are the same – one large Australian telco sells conference calling enabled handsets while their network does not support conference (e.g. 3 way) calling. This is not disclosed up front- we found out when we tried our first conf call. Ask the question, be wary of penguins and remember Skype is your friend.

Hope these few thoughts help. There is more we are learning each day so will stick up some more thoughts soon.

Cross-platform development

2010-10-04T11:36:00.003+11:00

During the course of developing Biarri’s flagship Workbench product, we’ve taken pains to ensure that our (GUI-less) optimisation “engines” work well under both Windows and Linux operating systems (so-called cross-platform). This turns out to be relatively easy as long as you stay away from the big OS-specific frameworks (e.g. Microsoft’s MFC/COM/ATL etc). We’ve picked up some handy tips along the way, particularly applicable to C++ development, which are worth sharing here.

Be aware of differences in line endings - Windows uses carriage return and line feed \r\n, while Linux/Unix uses just line feed \n. (Note that Visual Studio will show files with Linux line feeds correctly, but Notepad won't - this is one way to tell what line endings your file has in Windows). This can be particularly important when importing data e.g. into databases where the file originates from another OS.

Always use forward slashes for file paths, not backslashes. Also, file names and folder paths are case sensitive under Linux but not under Windows. And don't assume there is a C: or D: drive!

You may have to be careful writing to temporary files and folders. In Linux /tmp is often used; in Windows /[user]/AppData/local/temp (location of the TEMP environment variable – e.g. type “%TEMP%” into the start menu or Windows Explorer). For Linux, it is sometimes necessary to manipulate a folder’s “sticky bit” to ensure that the folder is accessible by other users (e.g. a Postgres database user) – e.g. in Python:

Be aware of the differences in file permissions in Windows and Linux. In Linux files have an “executable” bit. chmod a+x [file] makes a file an exe, which can then be run with “./filename”.

For C++ development:

Name all cpp and h files in lower case if possible. Files are case sensitive in Linux and this includes #include's!

For compiling with GCC under Linux, the last line in a C++ file must be blank.

In Linux C++ programs, general exception handling with catch(...) does not work. You can use sighandlers instead (see this for example), though it's not as good – it is more equivalent to an exit(), with a chance to clean up.

Beware doubles comparisons and inequality checking, at least in C++ programs. Always use a delta i.e. A == B may not be the case in both Windows and Linux if they are essentially the same number so use fabs(A - B) < delta, where delta is some suitably small number.

Build tips for Linux: Type "make" when you are in the directory to build the project. This will search for a file called "Makefile" and run it. (Use "make -f filename" to make from a different makefile). To force a recompile you can "touch" a file using "touch filename".
To clean out all object files type "make clean" (as long as your make file defines what cleaning does...). Use "make -j4" to run make with for concurrent jobs, to take advantage of multicore.

In bash, to get a recursive line count of .cpp/.h files: find [directory] -type f -name *.cpp -exec wc -l {} \; | awk '{total += $1} END{print total}'

Differences for GCC 4.3 compiler (e.g. compared to version 4.2.4): http://gcc.gnu.org/gcc-4.3/porting_to.html. (To check your version of GCC, type gcc -v).

Operations Research - Cloud Computing – The impending launch of Biarri’s WorkBench

2010-09-29T21:21:00.002+10:00

With the impending launch of Biarri’s workbench and our ongoing close relationship with Schweppes for the daily routing of soft drink deliveries (an application of perhaps the most well known operations research problem: the vehicle routing problem), I thought that the following excerpt from a journal article submitted to the Asia Pacific Journal of Operations Research would be a very timely blog post.

The journal article is entitled “Real-Life Vehicle Routing with Time Windows for Visual Attractiveness and Operational Robustness” and it describes the vehicle routing algorithm we have implemented for Schweppes.

The excerpt details a specific example encompassing two things we are very passionate about at Biarri. First “Commercial Mathematics” – that is making OR (well not strictly just OR) work in the real world. And second, the revolutionary capabilities that the advent of cloud computing has for the delivery of software.

“Vehicle routing problems manifest in a remarkably wide range of commercial and non-commercial enterprises. From: industrial waste collection to grocery delivery; underground mining crew replenishment to postal and courier collection and delivery; inbound manufacturing component transportation to finished car distribution; in-home primary health care delivery to pathology specimen clearances from surgeries for analysis; and from coal seam gas field equipment maintenance to beverage distribution, to name but a few.

Automated planning systems used by industry at present are predominantly client-server or desktop based applications. Such systems are often: expensive, requiring a large upfront capital investment; accompanied by a large software deployment project requiring initial and ongoing IT department cooperation; customisable to a particular organisations requirements, however commonly retain a large amount of exposed functionality due to the breadth of the existing client base; and require substantial user training as the workflow is usually not restricted in a linear fashion .... Each of these characteristics constitutes a barrier to adoption of automated planning systems, and for most small to medium enterprises these barriers prove insurmountable.

With the advent of cloud computing and software as a service (SaaS) these barriers are being removed. SaaS: embodies a different commercial model; has essentially no IT footprint; mandates (as vendors may never directly interact with potential clients) simple intuitive linear workflows; and involves almost no end user training beyond perhaps an optional demonstration video.

The emergence of this new avenue for the delivery of optimisation based planning systems heralds, a heretofore, unparalleled opportunity for operations research practitioners to engage with a wider potential consumer base than ever before. However, the nature of the delivery mechanism requires the algorithms developed: to be robust and flexible (within their domain of application they must be capable to dealing with a wide range of input data); to have very short run times (the user base is more likely to be under time pressure than ever before); to produce high quality solutions (noting the inherent trade off between run time and solution quality); to be wrapped in a simple linear workflow (meaning it is always obvious what the next step in the planning process is); but above all, be able to produce real-life, practically implementable solutions, without the need for user training and/or experience.

For pure delivery, or pure pick up vehicle routing applications, real-life, practically implementable solutions are often synonymous with geographically compact, non-overlapping routes with little or no intra-route cross over. There are numerous reasons why such solutions are preferred .... If a customer cannot be serviced at the preferred time (e.g. the vehicle cannot get access, the customer is closed, another delivery is taking place, the customer is too busy), because the route stays in the same geographical area, it is easy to return to the customer at a later time. During busy traffic periods drivers are loathe to exit and re-enter a motorway to service individual customers. Even though such customers may be enroute to the
bulk of the customers the route services, thus incurring a minimum of additional kilometres, they may nevertheless be far from the majority of the customers the route services. If there is severe traffic disruption, it is easier to use local alternate routes between customers in a route that is geographically compact to ensure that pick-ups or deliveries can still be made. Third party transport providers, which prefer routes to be as simple as possible, may exert some influence over the planning process. Finally ... it is easier to maintain customer relationships by assigning drivers to routes that routinely service a similar geographical area. In summary, solutions which are more visually attractive are more robust, and thus more likely to actually deliver the full extent of the cost savings that should flow from the use of automated planning systems.

This paper describes an algorithm for the vehicle routing problem with time windows, …. The algorithm is: robust and flexible; fast; wrapped in a user interface utilising a simple linear workflow and so requires no user training or experience; and produces high quality, visually attractive and practically implementable solutions.”

Our technology stack

2010-08-30T09:11:00.006+10:00

Over the course of developing our Workbench solution we've adopted a powerful set of interconnecting components. It’s worth mentioning what these are and how they fit together.

Almost all the components of the stack are free and/or open source. We want to be as platform independent as possible and not get too locked in to one technology paradigm. This means that as much as possible, parts should be as “hot swappable” as possible - which also helps encourage strong componentisation. Using components with mature and open/standardised interfaces is very necessary when you're crossing language boundaries (most notably, Javascript-Python and Python-C++) and client/server boundaries; otherwise you risk re-inventing the wheel. Ideally each component we use should also still be in active development (in the IT world - with the odd highly venerable exception - if software is not growing and evolving, it's usually either dying, already in it's death throes, or extinct).

There’s an art to using the right tool for the job, and we’ve made mistakes. We over-used Mako (see Loki’s blog post) and also originally used a slightly inferior lib for the C++ xmlrpc back end; both these mis-steps were fairly easily rectified. Arguably, we probably still use too much C++ and not enough Python – the C++ line count dwarfs the Python line count by a considerable margin. One last interesting point is that, at the moment, we're still eschewing use of an ORM (Object Relational Mapping layer - such as SQL Alchemy) - time will tell whether that is a good idea or not.

Server:

Python
- language

CherryPy
- Python-based web app framework

PostgreSQL
- open source RDBMS

PsycoPg2
- database adaptor for PostgreSQL/Python

XmlRPCLib
- XML RPC (used to communicate with some of the engines)

Mako
- Python template library

Repoze
- Zope/WSGI Python middleware (for authentication)

Client:

JavaScript
- client-side browser language

jQuery
- JavaScript library for event handling and more

Cloudmade
- OSM map provider/server

Data Interchange:

JSON - JavaScript Object Notation

XML - eXtensible Markup Language

Mathematical Engines:

Mostly in C++ using STL

CppUnit
- C++ library for unit testing

OGR
- map data library, part of GDAL - used to read map data

Libxmlrpc-c
- C++ back end for XMLRPC - used by a running process to communicate with the front end via Python

Logos

2010-08-09T14:59:00.006+10:00

The word that explains

Logos (λὀγοσ) is a Greek word that can be translated as "the word that explains”. Logo, as we have come to understand it means to produce meaning and which unifies a concept with a visual design. In ancient Greece, for Philosopher Heraclitus "Logos" was the place where the universal belongs. The universe can be seeing as chaos, as a river that flows constantly and never repeats the same pattern. But if we all head towards the "logos" we would achieve a universal way to understand the world. In order to see things in the same common way, we should "homo-logate" (ὁμολουεῖν) become one to the logos (homogeneous). 

In post-industrial times, Logo is the word that unifies a group of people working in a project under a name, usually represented through an icon, which is easy to identify. 
Heraclitus was worried about each one having its own logo so as to avoid chaos. In modern times, a business, a company, a project requires the same unification or "homogenisation" in order to ensure that all the people involve are on the same boat.  For a project to run properly, it needs "Logos".

When a graphic designer creates a logo, is not simply doing a nice drawing to promote a company; instead, it is discovering that universal symbol that belong to everyone involve in the project and unifies the project it self. There is nothing about aesthetics or make up; instead, it is the ultimate creative process, the beginning of the universal value that produces identification or "homogenisation" for the group. For this reason, an assembly line would never fit into an artistic conception; as in nature and life, there is no order in what should be first or second, instead, everything works and belong to a coherent collaborative process: The artist discovers the symbol that drives us all towards the universal and in this light it is where everyone can see themselves as belonging to the same thing.

  Another linguistic particularity is the word “project”, commonly used in our times. “Pro-eject” means to head towards something. If “logo” is the symbol that identifies the group, a “project” should be coherent with the logo, with its symbol. For example: an olive tree can identify a company that produces olive oil. There is a direct correlation in what this company does and its symbol. But we do not always find this obvious and direct relationship; usually, logos and promotional images are based on the concept of manipulation, distortion of the reality and many companies use their logo as an aesthetic camouflage. A very clear case of this is commonly identifiable in Petrol companies, most of them use green flowers or vivid environmental colours to mislead consumers into thinking their products are environmentally sound.

The idea of manipulation thought images evolved with the time responding to consumer’s expectations. Few decades ago, petrol companies where identified with symbols of power/energy while these days the new black is green.

Beyond moral, a designer is a communicator, not a manipulator. When Heraclitus invokes Greeks to follow the Logos (λὀγοσ) I guess he was hoping for a more creative world. A true artist should be a demiurge, someone capable of finding creativity between two worlds: the fluent chaotic river and universal beauty of the tree on the other side of the river.
That is what we become when we advocate to a pro-eject.

read more at: sebastianbourges.com

The 'L' Shaped Excel Graph

2010-08-06T10:44:00.007+10:00

I recently had to develop a graphical report in Excel that was required to fit perfectly on a single A4 page for printing purposes. A lot of information was supposed to fit on this page, so using the space well was very important. After the initial formatting was reviewed it was decided that one graph in particular needed to be bigger and clearer. With a bit of reshuffling the final space available for this graph was a big ‘L’ shape at the bottom of the page.

A perfect amount of space but unfortunately part of the ‘L’ would have to go to waste – unless some sort of ‘L’ shaped graph could be created…

The ‘L’ Shaped Graph:

The ‘L’ shaped graph utilises the extra space by moving the title and legend away from the actual data. This results in the data having space to expand and become more readable (very important with this detailed ‘Clustered-Stacked Column Graph’). It is important to note that once created, it is very difficult to edit this graph and it should most likely be treated as read only.

How it is done:

Finalise the data for your chart. (The legend will no longer update to reflect changes in the data once we are done)
Create a copy of the chart
Delete all axis labels and gridlines on the copy.

Now comes the trick – removing the data without affecting the legend:

To do this, select each data point on the copy (not data series) and right-click: Format data point. When you first click a data point it most likely selects the entire series, so click again to target just the data point. If you were to do this for the series then the legend would update to reflect the change, which we don’t want.
Select No-fill in the Fill tab. This will white out that data point. Do this for all the data points until just the legend and title remain on the copy.

Now to fix up the original chart:

Select the legend and title of the original chart and delete
Resize to utilise the space

The rest is just simple formatting to customise the graph to your needs. Here are some of the things I did:

Select the legend of the copy and right-click: Format legend
In Legend Options: Legend Position select ‘Top’.
Now resize the legend of the copied chart to whatever you like. You may notice that the chart area, although blank, is still selectable and may be an annoyance while resizing other aspects of the copy. This can not be deleted however you can resize it to be very tiny in one of the corners of that graph.
Remove the borders of both the original graph and the copy. Right-click: Format Chart area: Border Colour: No line.
Now create borders to surround the ‘L’ shaped chart on the spreadsheet itself to make it appear like the original and the copy are joined.

The process of whiting out each data point can be very tedious. In practice I created this graph with VBA to automate the process. There are some more tricks required if using VBA for the whole graph however I won’t go into them here. Like all VBA, heavy use of the macro recorder, trial and error, a little blind luck and a lot of patience go a long way.

An IDE for Python

2010-07-27T11:44:00.004+10:00

For some time now I have been trying to find a decent IDE with step-through debugging support for Python. I've wanted it for Linux, but Windows support would also be a bonus.

There’s some debate about the need for an IDE for Python, which (as a veteran of C++ development with Visual Studio) I am still pondering. I get that Python is a higher level language (umm, that’s why I’m using it), but the central problem of ironing out the kinks in the business/engine logic of my code is never going to go away. It really makes me wonder what size and types of code bases the IDE/debugger naysayers are building.

People also talk about Eclipse with PyDev, but I’m deterred by the reputedly formidable learning curve, the reportedly sluggish performance, and the apparent bloat of it. I wanted something lighter, but still free. And I didn’t want something that would require a big project hierarchy, settings tweakings etc, just to run a small Python script. I don’t think my needs are outlandish: easy thing should be easy, hard things should be possible...

This comparison of Python IDEs – the first hit on Google – seems good but is 5 years old (ancient in software development terms). And the Wikipedia comparison table is just the basic facts, ma’am. The Python.org list of IDEs is better, but without some sort of detail or commentary it's difficult to figure out what's best for your needs, and parts of it are out of date. Mile long feature lists are all well and good, but how well do the features do what they're supposed to do?

So I embarked on a trial of a few of the free IDEs out there. First stop was SPE – “Stani’s Python Editor”... which I couldn’t get installing. I know, I know, in Linux Land you’ve got to be prepared to tinker in the innards with spanner in hand... but a frustrating hour or two later, no go. Perhaps because this tool doesn’t seem to be actively developed any more (as I found out afterwards). Next I tried Boa Constructor. It installed, and first impressions were cautiously positive, though it felt like beta software to me. Sure enough after trying to use it in anger the pain points came – I couldn’t figure out the rhyme or reason why it wouldn’t just run the Python script I had open, I had to constantly restart, breakpoints weren’t always respected, etc. Overall it seems more aimed at GUI building than running scripts.

Next was the Eric IDE. Eric installs with just a simple “apt-get install eric”. The Python file you have open runs with “Start/Debug Script...” Breakpoints and stepping through just work (in fact, debugging is an absolute breeze). Lines with Python parse errors get a little bug icon on them in the editor margin (cute, but also handy). It’s not perfect by any means – it takes a while to start up, it occasionally automatically breaks at places where you have no breakpoints, etc etc. It’s GUI is formidable at first glance, but it’s set out logically and should seem familiar to those like me steeped in Visual Studio. It's also still being very actively developed.

One interesting aspect of Eric's editor is that it uses rich text with proportionally spaced fonts, in contrast to the plain monospaced font that most code editors sport. This might seem sacrilegious to some, but it seems to work fine for me, and in fact lets me see more code on the screen. It's not so good for ASCII art though.

Obviously I haven’t tried trials of the commercial IDEs out there – Komodo Edit, WingWare, etc, and I’d be curious what "killer features" (that work out of the box!) they have that Eric doesn’t. But for now, the journey’s over, with a clear winner.

Choosing a web framework and what to do about templates

2010-07-06T15:50:00.003+10:00

In the ongoing evolution of the workbench I've solved some interesting problems and found some hard ones I haven't solved yet. I also made some poor choices and later found more elegant solutions. All this is a learning process and I'd like to share some of those findings in this series. I'll cover templates, web framework choices and making your own widgets.

Templates:

Originally I searched for the most elegant and expressive templating language for what I needed to do and found Mako. It's fast, good for simple interpolation and some basic control structures. However, after trying to write even one of our basic workflows in it, I was less enthused. I needed a range of slightly dodgy hacks to get what I wanted out.

For instance, I wanted to have a mako function put a piece of text (javascript) in a buffer and output it later (in the head section). I found a way of doing this:



<%!
    def string_buffer(fn, target_buffer='html_buf'):
        def decorate(context, *args, **kw):
            context['attributes'][target_buffer]  = runtime.capture(context, fn, *args, **kw)
            return ''        
        return decorate
%>

<%def name="doc_ready_js()" decorator="partial(string_buffer, target_buffer='doc_ready_buf')">
    ${caller.body()}
</%def>
<%def name="main_js()" decorator="partial(string_buffer, target_buffer='main_js_buf')">
    ${caller.body()}
</%def>

Needing to use decorators, partial function application and weird context/global variable magic for quite simple features freaked me out enough and convinced me to hunt for a better solution. I wanted the simple to be easy and the hard to be possible. The template was making the the simple hard. I looked at a pile of templating languages and I still liked Mako the most out of them, but decided after reading this wonderful rant by by Tavis Rudd that the solution is obvious: Do simple stuff with templates if you need to, but mostly avoid them.

The clearest way to demonstrate the advantages of a pure python solution is with code, both sections of code generate the same simplified workflow:

The Mako:



<%namespace name="zui" file="zui.mako"/>
<%namespace name="importer" file="importer.mako"/>
<%namespace name="widgets" file="widgets.mako"/>

<%inherit file="workflow_wrapper.mako"/>
<%def name="workflow()">
    <%
        current_workflow_id = attributes['get_new_id']()
        attributes['workflow_steps'][current_workflow_id] = [[],[],[]]
        grid = attributes['get_new_id']()
        map = attributes['get_new_id']()
    %>

    <%call expr="importer.importer({'name':'text', 'address':'text'}, ['name', 'address'], r'''
    Here you need to select a CSV file with address and a name for the locations
    you want to geocode. Names have to be unique.''', table_name,  'import_done()')"></%call>

    ## events are run on the client side so are javascript
    <%call expr="widgets.run_engine_step('Geocode',
                                         {'onclick':'''calculate_button('geocode_addresses',
                                         geocode_results);'''},
                                         workflow_icon='geocode')"></%call>

    <%call expr="zui.workflow_step('output', 'View Results')">
        ${widgets.grid_holder(table_name, grid)}
        ${widgets.map_holder(map)}
    </%call>

    <%zui:main_js>
        function import_done(){
            console.log("run when import is finished");
        }

        function geocode_results(){
            ${widgets.grid_init(grid, table_name, '')}
        }
    </%zui:main_js>
</%def>

The Python:



import zui, widgets, importer, layout

class Geocoder(zui.Workflow):
    def __init__(self, the_zui, table):
        zui.Workflow.__init__(self, the_zui, table)
        importer.Importer(self, {'name':'text', 'address':'text'}, ['name', 'address'], r'''
            Here you need to select a CSV file with address and a name for the locations
            you want to geocode. Names have to be unique.''', self.table_name,
            'Geocoder.import_done()')

        ## events are run on the client side so are javascript
        widgets.RunEngineStep(self, 'Geocode',
                {'onclick':'''calculate_button('geocode_addresses',
                    Geocoder.geocode_results);'''}, workflow_icon='geocode')

        grid = widgets.Grid(self.table_name)
        cm_map = widgets.Map()

        self.workflow_step('output', 'View Results', container_size = "container-results1", body =
                layout.multi_column('', grid, cm_map))

        self.main_js(
            '''
            function import_done(){
                console.log("run when import is finished");
            }
            function geocode_results(){
                '''+grid.load()+'''
            }
            ''')

It is so much easier to be expressive and elegant when you minimise the usage of templates. You can clean things up by using some of the expressive power of a multi-paradigm language and not have to do anything too strange. I still see the use for templates when you need to interpolate a few variables into a long piece of HTML. The built in python template function seems fine for this. In my next post I'll over web frameworks and what I think about widgets. You should be able to see from today's bit of code that I think they are important and hard to do with templates.

Loki

p.s.

If you feel like relaxing after reading this, check out my band's YouTube channel. :)

Javascript: It's actually pretty decent.

2010-07-02T15:22:00.004+10:00

I admit it, the 'Java' in the word JavaScript has previously put me off. However, it turns out it really is rather nice and with some searching you can find some good documentation!

In our midst we have a language with more in common with scheme and lisp than Java or VB! And, it's wonderful features are eloquently explained by Douglas Crockford covering how to use encapsulation and a great series of slides by John Resig explaining closures and prototype based inheritance.

Douglas Crockford mentions a JavaScript gotcha that bit me - "this" which you use like self in python but is implicit, is set incorrectly for inner functions. It's simple to get around though, just set a var, i.e "that" to it's value before your inner functions.

I'm planning, and I admit this plan may be impeded by other more pressing concerns, to write a series of posts explaining my recent investigation of python (and some non-python) web frameworks. I also want to discuss my quest to write a application framework based on the principles of simplicity, elegance and power, for users and developers.

Despite JavaScript's positives the idea of a python to JS complier like Pyjamas does seem enticing. I'd be really interested in hearing comments from anyone who's used it. Python to me is still far more elegant. :)

Loki

Merging PostScript (.ps) files

2010-06-28T14:36:00.002+10:00

I recently had to merge the output of a bunch of academic papers I had written in Latex. Each file uses a document class from a range of academic journals (has its own title, abstract, bibliography etc).

I could have made a single Latex file and shoe horned each individual file into it, but I knew there had to be a better (less labour intensive) way.

Using Latex you can create pdf (using dvipdfm) or ps (using dvips) files. There are various ways to merge pdf files (pdfsam) and ps files (using ghostscript).

However the quality of the merged file I produced was always very poor. Eventually (after much googling) I found a solution

gswin32c.exe -dNOCACHE -dNOPAUSE -sDEVICE=pswrite -dBATCH -sOutputFile=Output.ps Input1.ps Input2.ps Input3.ps Input4.ps

The "-dNOCACHE" option preserves the quality of the output file. The output file however is very large.

To ensure the page numbers in the merged document are continuous you can use the Latex command "\setcounter{page}{X}"

Thoughts on Point Solutions

2010-06-28T11:33:00.002+10:00

Lately I have been thinking a bit about the advantages of small, tightly focussed web apps (so-called “point solutions”) that scratch a single little itch, versus larger, more powerful and general web apps that tend to deliver more of a total body rub. This question is of utmost importance to a company like Biarri that needs to place its development time and effort into the best channels.

The question was highlighted by a real-world problem a colleague posed recently: how to assign foursomes in rounds of Golf so that all of the players got to play with each other player at least once. It is not trivial to construct such a solution (if one even exists) by hand, if the constraints are “tight” enough (for example, 20 players and 8 rounds).

Small point solutions that solve a small but non-trivial problem like this might be fairly quick to develop and deploy on the web. But it doesn’t take much feature creep before you get a pile of extra “features” (particular requirements for some players, minimising the number of repeated pairings, right through to printing out score cards etc); before you know it (or more precisely, after months or years of hard coding) you’d have a full-blown Golf Tournament Scheduler. Such a web app might sell for much more, but would probably attract many less customers. And what happened to the poor casual golfer or golf tournament organiser on a shoestring budget who just wanted to solve his or her original golf player assignment problem?

In the spirit of acknowledging that the future is impossible to predict, I think Biarri must address more wide-ranging, lightweight “point solutions”, particularly at our fledgling stage. More mini-apps with a wider potential customer base will allow us to gauge which itches need the most scratching; more complex apps, as every seasoned developer knows, seem to always cause issues and problems – in short, sheer complexity – quite out of scale with the larger code line count; not to mention being harder to use and understand for users (more buttons!)

Those who have test-driven our Workbench solution will also know that, to some extent, we’re trying to have our cake and eat it to, by allowing these smaller “point” solutions to exist as workflows (standalone web apps) in their own right, whilst also being “nestable” – that is, able to be composited in a larger, more powerful workflow. Look out for Geocoding as a sub-workflow inside Travel Time Calculation, coming to the Workbench very soon. And who knows if the Biarri Golf Tournament Organiser will ever eventuate!

Re-solving Linear Programs with COIN-OR

2010-06-24T12:28:00.007+10:00

I recently had to go from solving a Linear Program (LP) once only, to modifying the problem and re-solving it multiple times. I was using the excellent open source COIN-OR libraries as my interface (OSI) and solver (CLP).

I found that you will get obscure-looking crashes (though, to be sure, I am using the downloaded binaries, not building COIN from source) inside OSI/COIN if you do not keep around the matrix and vector objects that you use to build your LP. It is safe to delete these objects if you are solving the problem once only, even before you have called initialSolve(). But it is not safe if you have to re-solve the LP!

A simple generic little helper class I used for this is:

class LPObjects
{
public:
CoinPackedMatrix* m_matrix;
double* m_objective, * m_col_lb, * m_col_ub, * m_rowrhs, * m_rowsen;

LPObjects()
{
m_objective = NULL;
   m_col_lb = m_col_ub = NULL;
   m_rowrhs = m_rowsen = NULL;
m_matrix = NULL;
}

void Create(int n_cols, int n_rows, double dblInfinity)
{
m_objective = new double[n_cols]; //the objective coefficients
m_col_lb = new double[n_cols]; //the column lower bounds
   m_col_ub = new double[n_cols]; //the column upper bounds
m_rowrhs = new double[n_rows]; //the row rhs
m_rowsen = new char[n_rows]; //the row sense

for (int i = 0; i
m_objective[i] = 0.0;
for (int i = 0; i
{
m_col_lb[i] = 0.0;
m_col_ub[i] = dblInfinity;
}
m_matrix = new CoinPackedMatrix(false, 0, 0);
m_matrix->setDimensions(0, n_cols);
}
};

If si is your OsiXxxSolverInterface object, and g_LP is a static pointer to an instance of LPObjects, then g_LP can be initialised with:

if (g_LP == NULL)
{
g_LP = new LPObjects;
g_LP->Create(n_cols, n_rows, si->getInfinity());
}

Create your problem representation by filling out rowrhs, rowsen, etc, and be sure to use CoinPackedMatrix’s appendRows to populate your matrix coefficients in bulk rather than using multiple calls to appendRow (which is much much slower for larger LPs).

Then load your problem into OSI with:

si->loadProblem(*g_LP->m_matrix, g_LP->m_col_lb, g_LP->m_col_ub, g_LP->m_objective, g_LP->m_rowsen, g_LP->m_rowrhs, NULL);

and solve with si->initialSolve(). Now you change the LP using the usual methods such as si->setColUpper(), and re-solve with si->resolve(). Note your call to si->setColUpper() (or whatever) does not directly reference the matrix or vectors of g_LP, but those objects must remain in memory!

Code used with binary distribution “CoinAll-1.2.0-win32-msvc9” under Windows Vista / Visual C++ Express 2008.

Animating the trajectory of a vehicle routing optimisation algorithm

2010-06-18T22:37:00.011+10:00

I recently wanted to view an animation of a vehicle routing optimisation algorithm I have been working on.

The algorithm uses a two phase optimisation approach. The first phase is a construction heuristic. The second phase is guided local search meta heuristic. The algorithms primary focus is to produce a visually attractive solution. Visually attractive solutions are usually synonymous with compact, non-overlapping routes with little or no intra-route cross over.

A colleague suggested I have a look at pygame. I hunted around to try and find an existing script that did something similar to what I wanted and found the fracture simulator. It is always easier to modify something than to start from scratch!

I was pleasantly surprised at how easy it turned out to be. I am very much a python novice, and had never heard of pygame before, but have been programming for a number of years. I was able to get the basics of my animation up and running in one evening! This far exceeded my expectations!

The python script can be viewed here

The input file used to create the uploaded animation can be viewed here

Open up a World

2010-06-04T11:18:00.014+10:00

This post is the start of a new series about design. There will be a new post every few days explaining the Biarri approach to design.  

  The role of designers in modern times.  

The idea of graphic design is usually associated with aesthetics, the make up that’s put on top of an already developed product or concept in order to improve its image and sales. Under this point of view, design is just a part of production, an extra adjustment in the assembly line to achieve the final goal: increase sales with a more seductive image. 

 A designer is a person with creative sensibility, often with an artistic background; and good skills in generating visual ideas. But the most important point: A designer is someone who understands the creative process from the beginning to bring things into being, into existence.

In the times of Leonardo da Vinci, it wouldn't be surprising to find an artist capable of painting canvas and conceive helicopters, or a sculptor who was at the same time an expert in anatomy. A designer needs to have a holistic approach as a creator, inspired by the discovery of the beauty and harmony present in every creative process.

Since the times of the industrial revolution the process of art has been displaced from being an independent conceptual process itself to become just a part of an assembly line. As industries grow, the role of the designer has become more and more specific. Currently we can find fashion designers, interior designers, industrial designers, graphic designers, communication designers, interface designers, light designers, set designers, etcetera; and within each category they can even become more specialised: there are designers specifically for typography, branding, packaging, magazines, print industries, concept artists, art directors, illustrators, 2D and 3D artists and so on.

The artist is still responsible for the idea generation but now from a specific position. Their work is no longer to bring something into existence; instead it is limited to bringing only the aesthetic into existence. The product or concept has already been determined, in many cases, by a business group that works closely with the financial and marketing department.  This changing role has serious implications in terms of the way we value the artistic process. The main issue is the frustration for the artist in loosing their creative autonomy, which is by nature a holistic process, ecstatic.

Vertigo and Ecstasy

Any creative person will be sensitive to experiences of “ecstasy”, from the Greek “existanai” which means: “drive out of one’s mind”, something that comes from inside to outside. Ecstasy is an experience that born to be share.
People consciously share ecstatic experiences during creative processes of any type. The result always opens up a new world of possibilities. Creativity then is a constant starting point to new possibilities.

On the other hand, we have experiences of “vertigo”. From Latin “vertere” originally used to represent the whirling movement of liquids. An experience of Vertigo opposes to the experience of ecstasy: “vertigo” is something that comes from outside to inside and is experimented individually. An experience of vertigo will not open up a world or any possibilities. It’s just a pleasant short moment, a feeling that ends up when the stimulus is not present or relevant anymore.

From the consumer point of view, an element of consume is an experience of vertigo. As its name indicates: consume, it extinguish it self once the short satisfaction is reached.
From the designer point of view, creativity applied to process of vertigo is destiny to die like most of the pages of a shallow magazine.

A real design should focus in open up possibilities for new creations. It is vital for all artist and people in general to create experiences of ecstasy in order to survive in time.

Ecstatic collaborative processes are the key to success.

Read more at: sebastianbourges.com

SQLite

2010-05-14T17:01:00.003+10:00

I recently needed to organise and filter a heap of data from a new client. I didn't want to deal with the overhead of a full-blown database and decided to try sqlite3. As it turns out, it was really easy to work with since the bindings are included with Python2.6. All I needed to do was read a bit on how to interface Python with sqlite here: http://docs.python.org/library/sqlite3.html#module-sqlite3 and I was good to go! I also installed a nice database management utility called SQLite Database Browser v2.0b 1 which you can get here: http://sqlitebrowser.sourceforge.net. It makes managing the structure of the database a bit easier than working in a windows command prompt and you can write SQL on-the-fly if you're having some problems with your Python. I find that it is pretty stable (though some of my poorly written SQL queries do send it into a tizzy and I need to kill it and reopen).

Now, I've decided to create a sqlite database and integrate it with the Excel front-end for one of our solvers (used when clients require desktop deployment). I anticipated that integration with an SQL database would greatly simplify and speed-up the reporting (with the added bonus of a significant reduction in the need for me to write complex vba code). Initially I banged around getting really frustrated with Excel, and DAO (even after I installed the ODBC driver available here: http://www.ch-werner.de/sqliteodbc/). Then I discovered SQLite for Excel here: http://sqliteforexcel.codeplex.com/. Whew! So far, I have found it very easy to work with and I am busy completing my reporting tool.

The Mighty Middle Click

2010-04-23T10:42:00.002+10:00

I recently discovered a nice set of shortcuts that use the middle mouse button click:

• In Firefox or Internet Explorer, middle-clicking on a link will open that link in a new tab.
• In Firefox or Internet Explorer, middle-clicking on a tab will close that tab. This also works to close files in Visual Studio.

I’ve started to use these and find myself using them all the time. Apparently middle-click is useful in Windows 7 as well.

For more juicy Firefox shortcuts:
• Ctrl-T for a new tab
• Ctrl-W to close the currently active tab
• Space and Shift-space to page up/down
• Ctrl-F to find
• Ctrl-Tab/Ctrl-Shift-Tab to go forward/back one tab

And my two best Windows shortcuts, which I can't live without:
• Ctrl-Shift-Esc – bring up Task Manager
• F2 in Explorer – rename a file

Thanks also to http://www.knowledgesutra.com/index.php/15-firefox-tricks_t49240.html for some nice tips.

XSLT vs Python/LXML

2010-03-23T12:04:00.002+11:00

I recently had to adapt some old XLST which read in an XML document and did some transformations to turn it into CSV data. For those who don’t know what XSLT is (consider yourself lucky!), it is a declarative, XML-based transformation language usually used for transforming a source XML into destination XML.

Now the XSLT I was dealing with also had some Javascript and VBScript functions inside it, and after struggling with it for a while I eventually realised – to my horror, as I needed it to run cross-platform on both Windows and Linux – that it also incorporated some Microsoft-specific extensions. So I ditched the XSLT and switched to writing it from scratch in Python with the LXML library instead. Less than 2 hours later, to do the exact same task – including thorough error checking – the Python turned out to be 205 lines, while the XSLT was 714 lines.

For a long time I wondered if it was just me, that I was too much of a procedural, C++ thinker, and just didn’t “get” XLST. XSLT is supposed to be a purpose-built tool for the job, right? Well, I'll say now what I've always secretly thought, namely that XSLT is obtuse and horrible and I’ll steer clear of it forever. And evidently I’m not alone.

Using Simpy and Python

2010-02-02T15:52:00.004+11:00

I recently completed a discrete-event simulation model using SimPy. This was my first foray into Python programming and the first time I used a non-graphical discrete-event simulation package (most of my previous experience was using Witness).

The model tracks the utilisation of wagons on trains. Key constraints include availability of staff, availability of a path through the network for movement of trains and wagons, and availability of a spare locomotive to make ancillary wagon movements. Maintenance activities are simulated as are wagon faults and repairs.

The rail network is represented using a .dot file; our brilliant computer scientist, Loki, devised a mechanism to find a path through the network using the pygraph minmax module (pydot is used to parse the .dot file).

Although the learning curve was quite steep for me, I found that I really enjoy working with both Python and SimPy. It was liberating not fiddling with graphics which must be placed precisely on the screen or manipulated to get element iteractions to work perfectly. The use of a .dot file to represent the rail network means that it can be changed at will without requiring a vast amount of re-work to fix graphics (the initial network is a high-level representation of the actual network which will be refined in the future). As an added bonus, Graphviz can be used to view and modify the .dot files. SimPy doesn't come loaded with pre-configured entities (like machines, staff, tracks, forklifts, etc.). However, you can create any entity you desire and define it's behaviour using the Python programming language. Since SimPy is open source all code is freely available and can be modified to provide further functionality. This provides far greater flexibility than I've experienced with other packages since I was not limited by behaviours defined by the package developers and the use of a package-specific macro language. As with other simulation packages I've used, I created variables to track important information and then created both .txt and .csv output files on which extensive analysis can be performed.

While the cost of using Python and SimPy was a significant learning curve (for me) I feel that it was definitely worth it. We'll be looking to use SimPy on future projects!

Special thanks to Klaus Muller for his patience in answering my questions and assistance!

Lightweight Geocoding

2010-01-28T11:24:00.004+11:00

At Biarri we're always on the look out for inexpensive, light tools that we can integrate with and pass on the value to our Workbench customers. We've been hunting for quite a while but haven't yet found a good free Geocoder that is suitable for possibly intensive commercial use (and preferably worldwide!). Geocoders turn addresses into latitude/longitude coordinates; reverse geocoding, which is also sometimes useful, finds the closest valid address given a coordinate.

Any geocoder worth its salt needs to be able to do a good job at fuzzy string matching because addresses are frequently given incorrectly (e.g. "Road" instead of "Street", multiple names for highways - "Hume Highway" vs "M31", common misspellings, ambiguity at suburb boundaries, confusion over locality names, etc). It's the quality of this matching in combination with the correctness and up-to-date-ness of the underlying map data that sets different geocoders apart.

I recently tried out TinyGeocoder. It has a simple API which appears to work quite well. It doesn't give you any direct information about what it matched to or how well it matched, though you can use the reverse geocoding function to get an idea. Though I'm finding the results a little bit patchy and fickle (some well known streets just inexplicably won't geocode, for example), and not particularly fast, overall it is a reasonable tool.

Generic Geometry Library

2010-01-27T17:25:00.004+11:00

To anyone who is trying to download the Generic Geometry Library (which was recently accepted by Boost) - you need to use the following subversion command line

svn co https://svn.boost.org/svn/boost/sandbox/ggl geometry

and not the one that is currently advertised

svn co https://svn.boost.org/svn/boost/sandbox/geometry geometry

We plan to use this library to compute convex hulls, and area of intersection of polygons, within our engines that are written in c++

Biarri Turns One

2010-01-05T13:27:00.000+11:00

Today is Biarri's Birthday! We started trading as a business on 5 January 2009. We have had a great first year and we would like to thank everyone that have supported us. We look forward to building on our success and delivering smart optimisation tools in 2010 and beyond.

Designing Distribution Networks for Utilities

2010-01-05T13:19:00.000+11:00

Many utilities face the problem of designing distribution networks. These networks need to reach all (or nearly all) households to distribute (or collect) the commodity supplied by the utility. This may be electricity, gas, water, data (cable or optical fibre) or sewerage.

Typically these distribution networks are built as a subset of some underlying network. Most commonly this is the street network. The reason for this is obvious – the street network connects all households and businesses – which are the final demand points of utilities. Some networks may piggyback on existing networks. For example an optical fibre rollout may be done using existing power poles.

From a mathematical point of view we can consider the underlying points of demand (households or small collections of households) as nodes in a graph. The potential connections between these nodes are arcs. We refer to this network of potential connections as the underlying network. If the distribution network can be represented as a tree (and this is often the case) then the problem of determining the best distribution network is similar to the minimum spanning tree problem solved on the underlying network (http://en.wikipedia.org/wiki/Minimum_spanning_tree). This problem is straightforward to solve with an algorithm that is relatively easy to understand. The main information required is the cost of the candidate arcs. For example, it may be relatively cheap to string new cable between existing poles (for electricity or data networks), but much more expensive to dig a trench to install cable underground in an area with no overhead poles.

The problem faced by a utility will however have many practical considerations that cannot be handled by the classical minimum spanning tree algorithm. These will include:

The problem is really one of building several spanning trees, each of which can handle a maximum demand. This arises because of underlying technical restrictions such as transformer capacity; the need to have pumping stations for water, gas and sewerage; the maximum number of optical fibre connections to a hub; etc.
The cost of arcs in the tree depends on the downstream demand. This arises because the pipes or cables need to get bigger as the tree collects more downstream demand.
Splitting nodes in the tree may be expensive. This can occur because of splicing costs (in optical fibre networks) or the need to distribute pressure evenly.
The spanning trees will themselves need to be connected back to the “mains” network. Indeed the problem of designing the mains network may need to be considered as part of the overall problem, though typically most of the cost will be in the feeder network due to the sheer number of connections required.
The maximum network length from the mains to a node is limited due to pressure/voltage losses.

Recently Biarri has undertaken a Proof of Concept project with a utility to look at designing a distribution network. We have formulated and solved an optimisation based approach to give solutions guaranteed to be within a very small margin of the optimal solution. On a sample network covering around 2000 households we were able to achieve:

An optimal solution in two minutes on standard PC’s
6% reduction in equipment costs
Huge reduction in planning effort – it took many man weeks to design the same area by hand

The speed of the solution process frees up the network designer to consider higher level issues that are not part of the optimisation model, such as the possible saving from using different technical approaches. These choices can be modelled as a series of “What if…?” scenarios.

We plan to extend the work done so far to cover even more of the practical issues faced by utilities. A technical paper covering the work will be prepared in 2010.

generating kml thematic maps with python and postgis

2009-12-29T12:30:00.004+11:00

The Biarri workbench now has a tool for showing the frequency of an event per postcode. The colourised map (funkier name Choropleth map) is generated as KML from postcode regions that are stored in a postgis database. This just requires a little bit of sql and python:



dict_cur.execute("select p.gid, p.poa_2006, ST_AsKML(the_geom) as \
                  polygon, frequency, (frequency - stats.min) / (stats.max - \
                  stats.min)::float as scaled_frequency from "+table+" c, \
                  postcode_regions p, (select max(frequency), min(frequency) \
                  from "+table+") as stats where p.poa_2006 = c.postcode")
return my_lookup.get_template("thematic_kml.mako").render(rows=dict_cur.fetchall())

The sql just normalises the column we want as our colour intensity and returns the KML for each postcode region.

The mako template:



<?xml version="1.0" encoding="utf-8"?>
<kml xmlns="http://earth.google.com/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom/">
    <Document>
    <name>Thematic Map</name>
        <description>generated with the Biarri workbench</description>
    % for row in rows:
    <Style id="${row['gid']}">
        <LineStyle>
          <color>40000000</color>
          <width>3</width>
        </LineStyle>
        <PolyStyle>
            <color>${intensity(row['scaled_frequency'])}</color>
            <fill>1</fill>
            <outline>1</outline>
        </PolyStyle>
      </Style>
      <Placemark>
        <name>${row['poa_2006']}</name>
        <description>
        ${row['frequency']}
        </description>
        <styleUrl>#${row['gid']}</styleUrl>
        ${row['polygon'].replace(" ", '\n')}
      </Placemark>
    % endfor
      </Document>
</kml>
<%!
import webcolors
import colorsys

def intensity(frequency):
#frequencies are already scaled between 0 and 1
    return 'cc' webcolors.rgb_to_hex(tuple((int(a * 256) for a in colorsys.hsv_to_rgb(0.134, frequency, 0.90)))).strip('#')
%>

The mako template is quite easy to read if you notice the % sections are actually python code. The only challenge I had was that ST_AsKML produces kml with the co-ordinates space separated. Google earth is fine with this but the kml viewer we're using, web map lite, only likes them newline separated. My colour function just changes the intensity of the colour. A gradient between blue and red would probably be more appropriate.

Loki

Biarri's first Open Source contribution

2009-12-08T15:53:00.002+11:00

Today Biarri open sourced its first code. We've used Google Code hosting with an MIT License to open source a small program called osm2mif that converts Open Street Map data (OSM) to MapInfo Interchange Format (MID/MIF).

From the Google Code project homepage for osm2mif:

Converts an Open Street Map file (.OSM) to mid/mif output, creating a topology for routing (that is, breaking up the OSM ways at intersection points, and also processing "no right turn" restriction relations).

Can use filtering to convert only within a given coordinate bounding box; also can convert only a subset of features (e.g. just roads, or just coastlines and lakes/rivers, or some given mixture); can also transform tag values to given strings, and use particular MIF styles (pen style, polyline vs region etc) based on the given tag values.

The program is written in C++ and contains only one file, with less than 1000 lines, and no external dependencies apart from the standard library. Compiles under both Windows (Visual Studio) and Linux (GCC). Runs pretty quickly for OSM files under a Gigabyte, but does use up lots of RAM when running.

Sebastian Joins The Team

2009-11-29T15:41:00.000+11:00

We welcome Sebastian Bourges to Biarri as our new Artistic Director. His focus will be the design of the website and the usability of our optimisation tools. He will drive the online look and feel of the business.

Sebsatian is originally from Argentina where he studied Industrial Design at the National University of Cordoba. He has had a range of design roles including running his own graphic and multimedia design business. He migrated to Australia in 2006 and worked for Katana Productions which later merged with The Feds, where he developed a variety of multimedia projects and websites for film directors. Recently he has been working as an independent multimedia designer.

Sebastian is married, lives in Melbourne and is a keen photographer.

Using Open Street Map data for Australia

2009-11-27T23:01:00.004+11:00

I’ve been using Open Street Map (OSM) data for Australia for generating matrices of travel times and distances that are suitable for input to truck routing optimisation problems. OSM is a world-wide map data set that anyone can edit and is impressively comprehensive.

Some of the resources I’ve used are www.osmaustralia.org, which does a regular batch job to extract country by country OSM data, and http://keepright.x10hosting.com/ which keeps an up to date list of map data errors. These types of errors are important to know about when producing routable data from OSM, in particular dead end one way streets and “almost junctions”.

I’ve been using Quantum GIS (QGIS) to look at the data after converting it to MID/MIF format. I can get it to label the lines that are one way streets, but unfortunately there is no way to show the line direction in QGIS, which is a major pain. (Actually there is a way but you have to get a source code branch from their repo and build it all etc etc). [Now that I think of it, a silly cheat way to do it would also be to produce another field in the map data, say a letter which represents whether the arc is "D" for down or "L" for left etc (based on the difference in the lat/longs of the start and end nodes of each arc), and label that way.]

I also found that the performance of QGIS is quite different depending on whether the data is opened as TAB or MIDMIF format. TAB format is fairly fast (just zooming, panning etc) but MIDMIF is quite noticeably slower. You'd think it wouldn't make a difference as it would be using the same internal data representation, but obviously not.

I’m using some extra layers to show some of the processed data (picture below). For example, I have a layer which just shows the arcs involved in restricted turns, which I can layer on top of the street network and use a thicker line style for. I also use dots in another layer which I have produced based on my "island" processing. The dots represent "orphaned" nodes which are not on the main "island" of map data (which is connected in the sense that all nodes can reach all other nodes). These orphaned nodes will be ignored by the travel time calculation. There are around 9600 of these nodes in the entire set of Australia OSM routable data I’m using (which has 620042 nodes and 1503668 arcs in total). This filtered subset of OSM data I am using includes only street segments with a “highway” tag – this will exclude cycleways, pathways, hiking trails, ferry routes, ski slopes, building outlines, administrative boundaries, waterways etc.

Some observations on the OSM data (for Australia):
• One way information seems quite complete.
• Only 12% of the streets have road speed information (“max_speed” tag). This is an issue as vehicle routing needs highways to have a faster speed otherwise they won't be used (and there will be lots of rat-running for example). In longer routes (e.g. interstate) the travel time will also be grossly overestimated. A couple of things we could do here: search for all segments with "highway", "motorway", or "freeway" in their names and assume some sort of speed like 80km/h on these segments. Or, with a bit more coding effort, if there are chains of segments where one has a speed and the others with the same street name don't, assume that speed on all of those segments.
• About 70% of the road segments have street names. Some of the unnamed segments are bits of roundabout, service roads, motorway on/off-ramps etc. But, there are also a lot of streets which are classified (by their "highway" tag) as "secondary", "tertiary", "unclassified" and even "residential" which are not named. This is an issue when vehicle routing needs to produce driver directions or verbose path information.
• There's only several hundred instances of streets which have restricted manoeuvre information (banned right hand turns being the prime example). Most of these instances are in Sydney. In reality this number should be in the thousands or tens of thousands. This will likely be the biggest issue from a routing quality point of view - it will mean you can do many illegal right turns.
• I found a weird character or two at the end of the file, which causes both Python and C++ file reading functions to get confused and read blanks forever. Weird, but easily avoided.

Drawing figures for scientific publications: XFig and MS Windows

2009-11-18T11:15:00.002+11:00

I recently had to modify a picture (figure) as part of a resubmission to an academic journal.

The picture was created using XFig some time ago. XFig is one of the most widely used applications for creating figures for academic publications, and easily the best application I have used for this purpose. Its ability to incorporate Latex scientific formulas (use "xfig -specialtext -latexfonts -startlatexFont default" when launching XFig) and fonts, means that figures can be easily and beautifully integrated into Latex documents.

As I now have a windows machine, and am definitely no expert on linux (I am slowly trying to remedy this) I was dreading having to jump through the hoops to get XFig working on my new windows box (See http://www.cs.usask.ca/~wew036/latex/xfig.html). Part of the process involves installing Cygwin (a Linux-like environment for Windows). Not being familiar with linux this process seems quite convoluted (and I have been down this path before).

Googling XFig also brings up WinFig. Which is supposed to be very similar to XFig but runs on MS Windows. After downloading WinFig I quickly found out that you can only save figures with 15 or less objects in them without paying for the full version (making the free version not very useful). Something that the homepage neglects to mention.

I soon realised that because I had already installed virtual box and installed ubuntu (the process was very pain free) I should definitely try to use XFig within Ubuntu. Installing XFig within Ubuntu is what one my colleagues would call automagical - with a terminal inside Ubuntu type "sudo apt-get install xfig", then as mentioned before type "xfig -specialtext -latexfonts -startlatexFont default" and I was cooking straight away.

In order to open my file in XFig in Ubuntu I was hoping to be able to share some folders with MSWindows, and then mount them inside Ubuntu. Alas, despite all my efforts I have still not been able to get this to work. Email to the rescue, emailed them to myself, opened my gmail in ubuntu, saved the file and the problem was solved.

The moral of the story (well this blog) is that if you are trying to get XFig working on windows, - don't. Use the power that virtual box gives you and run XFig in Ubuntu within virtual box on your Windows machine! Now to move back to the Latex editing applications in linux and away from those I have been using with MS Windows!

Well Done Nick and Daryl - Maths Prize Winners

2009-11-17T14:00:00.001+11:00

Daryl Bruce and Nick Vaskrsic both won awards at the recent RMIT Mathematical Sciences Awards Night. Daryl won the Operations Research Prize and Nick won the Mathematics Rising Star Prize. Both of them currently work as part-time business analysts and route optimisation planners for Biarri.

It is great to see the guys doing so well in their studies. All the best to everyone in the Biarri team who are currently in the midst of exams - good luck!

Web Applications Designer

2009-11-07T09:21:00.002+11:00

We are making great progress with the Smart Optimisation Tools we are building for the Workbench. But we now need to pump up our web and application design capability so we can speed up development and polish the products. We are looking to hire a web application designer.

The job is based in Melbourne and we need someone who can design the 'look and feel' of our web applications, website and brand. This is a chance to get someone in on the ground floor to be the creative force behind our online style and usability. They will be Biarri's design guru. As a creative designer, they will have:

A strong visual eye and understanding of web usability
Creative flair that they use to build clean and appealing interfaces
The ability to drive the design process
End-to-end design skills, from sketches and wireframes to HTML/CSS integration

The job is also on Seek http://www.seek.com.au/job/industrial-web-applications-designer/melbourne/16207367/23/1/