October 10, 2015

Vince Knight

Getting your first lectureship

On Wednesday I was invited to participate in a Webinar entitled “Getting your first lecturship”. This was organised by Cardiff University’s graduate college (UGC). The format included an overview of recent findings from a survey conducted by the Association of Careers Advisory Services and then a discussion including Dr Sophie Coulombeau and Dr Claire Shaw.

You can see the recording below:

Here is the UGC page which includes a pdf with a variety of hopefully useful resources.

Sophie and Claire had some excellent advice and I would really recommend watching the video if you are at a stage of your academic career (PhD, postdoc…) and thinking of getting a lectureship in the United Kingdom (it’s probably useful in other countries also, but we concentrated on particularities of the UK).

Some of the particular points that stayed with me:

  • Self care: it’s very easy to let academia do a lot of damage. This is not something I’m particularly good at myself but stress and mental health is something that you should always keep an eye on: at whatever stage of your career you are.
  • Selling the whole narrative: at interview it’s important to find a way to think about what a particular post is looking for and how/why you fit that particular role.

What was also interesting was that we all seemed to agree that we were lucky. This was certainly the case with me, circumstances were just right and I was very lucky to have leadership that gave me a wealth of opportunities.

On a slightly related note, as I was writing this post, I thought about Tim Hopper’s great resource shouldigetaphd.com/. If you are thinking about doing a PhD at all, perhaps take a look at the awesome case stories there…

October 10, 2015 12:00 AM

September 30, 2015

William Stein

What is SageMath's strategy?

Here is SageMath's strategy, or at least what my strategy toward SageMath has been for the last 5 years.

Diagnose the problem

Statement of problem: SageMath is not growing.


Facts: Growth in the number of active users [1] of SageMath has stalled since about 2011 (as defined by Google analytics on sagemath.org). From 2008 to 2011, year-on-year growth was about 50%, which isn't great. However, from 2011 to now, year-on-year growth is slightly less than 0%. It was maybe -10% from 2013 to 2014. Incidentally, number of monthly active users of sagemath.org is about 68,652 right now, but the raw number isn't as import as the year-to-year rate of change.

I set an overall mission statement for the Sage project at the outset, which was is to be a viable alternative to Magma, Maple, Mathematica and Matlab. Being a "viable alternative" is something that holds or doesn't for specific people. A useful measure of this mission then is whether or not people use Sage. This is a different metric than trying to argue from "first principles" by making a list of features of each system, comparing benchmarks, etc.

Guiding policies

Statement of policy: focus on undergraduate students in STEM courses (science, tech, engineering, math)


In order for Sage to start growing again, identify groups of people that are not using Sage. Then decide, for each of these groups, who might find value in using Sage, especially if we are able to put work into making it easier for them to benefit from Sage. This is something to re-evaluate periodically. In itself, this is very generic -- it's what any software project that wishes to grow should do. The interesting part is the details.
Some big groups of potential future users of Sage, who use Sage very little now, include
  • employees/engineers in various industries (from defense contractors, to finance, to health care to "data science").
  • researchers in area of mathematics where Sage is currently not popular
  • undergraduate students in STEM courses (science, tech, engineering, math)
I think by far the most promising group is "undergraduate students in STEM courses". In many cases they use no software at all or are unhappy with what they do use. They are extremely cost sensitive. Open source provides a unique advantage in education because it is less expensive than closed source software, and having access to source code is something that instructors consider valuable as part of the learning experience. Also, state of the art performance, which often requires enormous dedicated for-pay work, is frequently not a requirement.


  • (a) Make access to Sage as easy as possible.
  • (b) Encourage the creation of educational resources (books, tutorials, etc.) that make using Sage for particular courses as easy as possible.
  • (c) Implement missing functionality in Sage that is needed in support of undergraduate teaching.


Why don't more undergraduates use Sage? For the most part, students use what they are told to use by their instructors. So why don't instructors chose to use Sage? (a) Sage is not trivial to install (in fact it is incredibly hard to install), (b) There are limited resources (books, tutorials, course materials, etc.) for making using Sage really easy, (c) Sage is missing key functionality needed in support undergraduate teaching.

Regarding (c), in 2008 Sage was utterly useless for most STEM courses. However, over the years things changed for the better, due to the hard work of Rob Beezer, Karl Dieter, Burcin Erocal, and many others. Also, for quite a bit of STEM work, the numerical Python ecosystem (and/or R) provides much of what is needed, and both have evolved enormously in recent years. They are all usable from Sage, and making such use easier should be an extremely high priority. Related -- Bill Hart wrote "I recently sat down with some serious developers and we discussed symbolics in Sage (which I know nothing about). They argued that Sage is not a viable contender in that area, and we discussed some of the possible reasons for that. " The reason is that the symbolic functionality in Sage is motivated by making Sage useful for undergraduate teaching; it has nothing to do with what serious developers in symbolics would care about.

Regarding (b), an NSF (called "UTMOST") helped in this direction... Also, Gregory Bard wrote "Sage for undergraduates", which is exactly the sort of thing we should be very strongly encouraging. This is a book that is published by the AMS and is also freely available. And it squarely addresses exactly this audience. Similarly, the French book that Paul Zimmerman edited is fantastic for France. Let's make an order of magnitude similar resources along these lines! Let's make vastly more tutorials and reference manuals that are "for undergraduates".

Regarding (a), in my opinion the most viable option that fits with current trends in software is a full web application that provides access to Sage. SageMathCloud is what I've been doing in this direction, and it's been growing since 2013 at over 100% year on year, and much is in place so that it could scale up to more users. It still has a huge way to go regarding user friendliness, and it is still losing money every month. But it is a concrete action toward which nontrivial effort has been invested, and it has the potential to solve problem (a) for a large number of potential STEM users. College students very often have extremely good bandwidth coupled with cheap weak laptops, so a web application is the natural solution for them.

Though much has been done to make Sage easier to install on individual computers, it's exactly the sort of problem that money could help solve, but for which we have little money. I'm optimistic that OpenDreamKit will do something in this direction.

[I've made this post motivated by the discussion in this thread.  Also, I used the framework from this book.]

by William Stein (noreply@blogger.com) at September 30, 2015 09:56 PM

September 22, 2015

William Stein

SageMathCloud's poor user retention rate

Poor retention rate

Many people try SageMathCloud, but only a small percentage stick around.  I definitely don't know why. Recent SageMathCloud rates are below 4%:

Is it performance?

Question: Are the people who try SMC discouraged by performance issues?

I think it's unlikely many users are leaving due to hitting noticeable performance issues.  I think I would know, since there's a huge bold messages all over the site that say "Email help@sagemath.com: in case of problems, do not hesitate to immediately email us. We want to know if anything is broken!"    In the past when there have been performance or availability issues -- which of course do happen sometimes due to bugs or whatever -- I quickly get a lot of emails.  I haven't got anything that mentioned performance recently.  And usage of SMC is at an all time high: in the last day there were 676 projects created and 3500 projects modified -- which is significantly higher than ever before since the site started.  It's also about 2.2x what we had exactly a year ago.    

Is it the user interface?

Question: Is the SMC user interface highly discouraging and difficult to use?

My current best guess is that the main reason for attrition of our users is that 
they do not understand how to actually use SageMathCloud (SMC), and the interface doesn't help at all.   I think a large number of users get massively confused and lost when trying to use SMC.    It's pretty obvious this happens if you just watch what they do...    In order to have a foundation on which to fix that, the plan I came up with in May was to at least fix the frontend implementation so that it would be  much easier to do development with -- by switching from a confusing  mess of jQuery soup, e.g., 2012-style single page app development -- to Facebook's new React.js approach.  This is basically half done and deployed, and I'm going to work very hard for a while to finish it.   Once it's done, it's going to be much easier to improve the UI to  make it more user friendly.

Is it the open source software?

Question: Is open source mathematical software not sufficiently user friendly?

Fixing the UI probably won't help so much with improving the underlying open 
source mathematical software to be friendly though.    This is a massive, deep, and very difficult problem, and might be why growth of Sage stopped in 2011:

SageMath (and maybe Numpy/Scipy/IPython/etc.) are not as user friendly as Mathematica/Matlab.   I think they could be even more user friendly, but it's highly unlikely as long as the developers are mostly working on SageMath in their spare time as part of advanced research projects (which have little to do with user friendliness).  

Analyzing data about mistakes, frustation, and issues people actually have with  real worksheets and notebooks could also help a lot with directing our  effort in improving Sage/Python/Numpy/etc to be more user friendly.

Is it support?

Question: Are users frustrated by lack of interactive support?

Having integrated high-quality support for users inside SMC, in which we help 
them write code, answer questions, etc., could help with retention.  

Why don't you use SageMathCloud?

I've been watching this  stuff closely for over a decade most waking moments, and everybody likes to complain to me.    Why don't you use SageMathCloud?   Tell me:  wstein@sagemath.com.

by William Stein (noreply@blogger.com) at September 22, 2015 07:52 PM

September 21, 2015

Sébastien Labbé

There are 13.366.431.646 solutions to the Quantumino game

Some years ago, I wrote code in Sage to solve the Quantumino puzzle. I also used it to make a one-minute video illustrating the Dancing links algorithm which I am proud to say it is now part of the Dancing links wikipedia page.


Let me recall that the goal of the Quantumino puzzle is to fill a \(2\times 5\times 8\) box with 16 out of 17 three-dimensional pentaminos. After writing the sage code to solve the puzzle, one question was left: how many solutions are there? Is the official website realist or very prudent when they say that there are over 10.000 potential solutions? Can it be computed in hours? days? months? years? ... or light-years!? The only thing I knew was that the following computation (letting the 0-th pentamino aside) never finished on my machine:

sage: from sage.games.quantumino import QuantuminoSolver
sage: QuantuminoSolver(0).number_of_solutions()   # long time :)

Since I spent already too much time on this side-project, I decided in 2012 to stop investing any more time on it and to really focus on finishing writing my thesis.

So before I finish writing my thesis, I knew that the computation was not going to take a light-year, since I was able to finish the computation of the number of solutions when the 0-th pentamino is put aside and when one pentamino is pre-positioned somewhere in the box. That computation completed in 4 hours on my old laptop and gave about 5 millions solutions. There are 17 choices of pentatminos to put aside, there are 360 distinct positions of that pentamino in the box, so I estimated the number of solution to be something like \(17\times 360\times 5000000 = 30 \times 10^9\). Most importantly, I estimated the computation to take \(17\times 360\times 4= 24480\) hours or 1020 days. Therefore, I knew I could not do it on my laptop.

But last year, I received an email from the designer of the Quantumino puzzle:

-------- Message transféré --------
Sujet : quantumino
Date : Tue, 09 Dec 2014 13:22:30 +0100
De : Nicolaas Neuwahl
Pour : Sebastien Labbe

hi sébastien labbé,

i'm the designer of the quantumino puzzle.
i'm not a mathematician, i'm an architect. i like mathematics.
i'm quite impressed to see the sage work on quantumino, also i have not the
knowledge for full understanding.

i have a question for you - can you tell me HOW MANY different quantumino-
solutions exist?

ty and bye

nicolaas neuwahl

This summer was a good timing to launch the computation on my beautiful Intel® Core™ i5-4590 CPU @ 3.30GHz × 4 at Université de Liège. First, I improved the Sage code to allow a parallel computation of number of solutions in the dancing links code (#18987, merged in a Sage 6.9.beta6). Secondly, we may remark that each tiling of the \(2\times 5\times 8\) box can be rotated in order to find 3 other solutions. It is possible to gain a factor 4 by avoiding to count 4 times the same solution up to rotations (#19107, still needs work from myself). Thanks to Vincent Delecroix for doing the review on both ticket. Dividing the estimated 1024 days of computation needed by a factor \(4\times 4=16\) gives an approximation of 64 days to complete the computation. Two months, just enough to be tractable!

With those two tickets (some previous version to be honest) on top of sage-6.8, I started the computation on August 4th and the computation finished last week on September 18th for a total of 45 days. The computation was stopped only once on September 8th (I forgot to close firefox and thunderbird that night...).

The number of solutions and computation time for each pentamino put aside together with the first solution found is shown in the table below. We remark that some values are equal when the aside pentaminoes are miror images (why!?:).

/Files/2015/b0.png /Files/2015/b1.png
634 900 493 solutions 634 900 493 solutions
2 days, 6:22:44.883358 2 days, 6:19:08.945691
/Files/2015/b2.png /Files/2015/b3.png
509 560 697 solutions 509 560 697 solutions
2 days, 0:01:36.844612 2 days, 0:41:59.447773
/Files/2015/b4.png /Files/2015/b5.png
628 384 422 solutions 628 384 422 solutions
2 days, 7:52:31.459247 2 days, 8:44:49.465672
/Files/2015/b6.png /Files/2015/b7.png
1 212 362 145 solutions 1 212 362 145 solutions
3 days, 17:25:00.346627 3 days, 19:10:02.353063
/Files/2015/b8.png /Files/2015/b9.png
197 325 298 solutions 556 534 800 solutions
22:51:54.439932 1 day, 19:05:23.908326
/Files/2015/b10.png /Files/2015/b11.png
664 820 756 solutions 468 206 736 solutions
2 days, 8:48:54.767662 1 day, 20:14:56.014557
/Files/2015/b12.png /Files/2015/b13.png
1 385 955 043 solutions 1 385 955 043 solutions
4 days, 1:40:30.270929 4 days, 4:44:05.399367
/Files/2015/b14.png /Files/2015/b15.png
694 998 374 solutions 694 998 374 solutions
2 days, 11:44:29.631 2 days, 6:01:57.946708
1 347 221 708 solutions  
3 days, 21:51:29.043459  

Therefore the total number of solutions up to rotations is 13 366 431 646 which is indeed more than 10000:)

sage: L = [634900493, 634900493, 509560697, 509560697, 628384422,
628384422, 1212362145, 1212362145, 197325298, 556534800, 664820756,
468206736, 1385955043, 1385955043, 694998374, 694998374, 1347221708]
sage: sum(L)
sage: factor(_)
2 * 23 * 271 * 1072231
The machine (4 cores) Intel® Core™ i5-4590 CPU @ 3.30GHz × 4 (Université de Liège)
Computation Time 45 days, (Aug 4th -- Sep 18th, 2015)
Number of solutions (up to rotations) 13 366 431 646
Number of solutions / cpu / second 859

My code will be available on github.

About the video on wikipedia.

I must say that the video is not perfect. On wikipedia, the file talk page of the video says that the Jerky camera movement is distracting. That is because I managed to make the video out of images created by .show(viewer='tachyon') which changes the coordinate system, hardcodes a lot of parameters, zoom properly, simplifies stuff to make sure the user don't see just a blank image. But, for making a movie, we need access to more parameters especially the placement of the camera (to avoid the jerky movement). I know that Tachyon allows all of that. It is still a project that I have to create a more versatile Graphics3D -> Tachyon conversion allowing to construct nice videos of evolving mathematical objects. That's another story.

by Sébastien Labbé at September 21, 2015 02:55 PM

September 17, 2015

Vince Knight

A Summer of game theory software development

This Summer has seen 3 undergraduates carry out 8 week placements with me developing further game theoretic code in Sagemath:

Hannah Lorrimore (going in to her 2nd year) spent her placement working very hard to implement classes for extensive form games. These are mainly a graphical representation of games based on a tree like structure. Hannah not only developed the appropriate data structures but also worked hard to make sure the graphics looked good. You can see an example of the output below:

James Campbell (going in to an industrial placement year) picked up where he left off last Summer (James built the first parts of Game Theory code for Sagemath) and developed a test for degeneracy of games. This involves building a corresponding linear system for the game and testing a particular condition. James and I wrote a blog post about some of the theory here: http://vknight.org/unpeudemath/code/2015/06/25/on_testing_degeneracy_of_games/.

Rhys Ward (going in to his first year) has been working at the interface between extensive form games and normal form games. His main contribution (Rhys is still working as of the writing of this) has been to build code that converts an extensive form game to a normal form game. This requires carefully traversing the underlying tree and keeping track of the strategy space. Rhys has also built a catalog of normal form games and is now starting to work on the capability to remove dominated strategies from a normal form game.

Hannah, Rhys and James have also been working in conjunction with Tobenna Peter Igwe who is a PhD student at the University of Liverpool. Tobenna has been implementing a variety of game theoretic code as part of the Google Summer of Code project with me as his mentor.

Hannah, James, Tobenna and I visited Oxford University to spend two days working with Dr Dima Pasechnik and giving a talk. You can see a video of the talk here: https://www.youtube.com/watch?v=v4kKYr5I2io

All of this code will now be reviewed by the Sagemath community and will (just as James’s code last year) be eventually available to anyone who wants to study game theory.

Note: this blog post is based on a similar Cardiff University newsletter item.

September 17, 2015 12:00 AM

September 10, 2015

William Stein

Funding Open Source Mathematical Software in the United States

I do not know how to get funding for open source mathematical software in the United States. However, I'm trying.

Why: Because Sage is Hobbling Along

Despite what we might think in our Sage-developer bubble, Sage is hobbling along, and without an infusion of financial support very soon, I think the project is going to fail in the next few years. I have access to Google analytics data for sagemath.org since 2007, and there has been no growth  in active users of the website since 2011:

Something that is Missing

The worse part of all for me, after ten years, is seeing things like this email today from John Palmieri, where he talks about writing slow but interesting algebraic topology code, and needing help from somebody who knows Cython to actually make his code fast.

I know from my three visits to the Magma group in Sydney that such assistance is precisely what having real financial support can provide. Such money makes it possible to have fulltime people who know the tools and how to optimize them well, and they work on this sort of speedup and integration -- this "devil is in the details" work -- for each major contribution (they are sort of like a highly skilled version of a journal copy editor and referee all in one). Doing this makes a massive difference, but also costs on the order of $1 million / year to have any real impact. 1 million is probably the Magma budget to support around 10 people and periodic visitors, and of course like 1% of the budget of Matlab/Mathematica. Magma has this support partly because Magma is closed source, and maintains tight control on who may use it.

Searching for a Funding Model

Sage is open source and freely available to all, so it is of potential huge value to the community by being owned by everybody and changeable. However, those who fund Magma (either directly or indirectly) haven't funded Sage at the same level for some reason. I can't make Sage closed source and copy that very successful funding model. I've tried everything I can think of given the time and resources I have, and the only model left that seems able to support open source is having a company that does something else well and makes money, then using some of the profit to fund open source (Intel is the biggest contributor to Linux).

SageMath, Inc.

Since I failed to find any companies that passionately care about Sage like Intel/Google/RedHat/etc. care about Linux, I started one. I've been working on SageMathCloud extremely hard for over 3 years now, with the hopes that at least it could be a way to fund Sage development.

by William Stein (noreply@blogger.com) at September 10, 2015 07:38 PM

September 05, 2015

William Stein

The Simons Foundation and Open Source Software

Jim Simons

Jim Simons is a mathematician who left academia to start a hedge fund that beat the stock market. He contributes back to the mathematical community through the Simons Foundation, which provides an enormous amount of support to mathematicians and physicists, and has many outreach programs.

SageMath is a large software package for mathematics that I started in 2005 with the goal of creating a free open source viable alternative to Magma, Mathematica, Maple, and Matlab. People frequently tell me I should approach the Simons Foundation for funding to support Sage. For example:
Jim Simons, after retiring from Renaissance Technologies with a cool 15 billion, has spent the last 10 years giving grants to people in the pure sciences. He's a true academic at heart [...] Anyways, he's very fond of academics and gives MacArthur-esque grants, especially to people who want to change the way mathematics is taught. Approach his fund. I'm 100% sure he'll give you a grant on the spot.

The National Science Foundation

Last month the http://sagemath.org website had 45,114 monthly active users. However, as far as I know, there is no NSF funding for Sage in the United States right now, and development is mostly done on a shoestring in spare time. We have recently failed to get several NSF grants for Sage, despite there being Sage-related grants in the past from NSF. I know that funding is random, and I will keep trying. I have two proposals for Sage funding submitted to NSF right now.

Several million dollars per year

I was incredibly excited in 2012 when David Eisenbud invited me to a meeting at the Simons Foundation headquarters in New York City with the following official description of their goals:
The purpose of this round table is to investigate what sorts of support would facilitate the development, deployment and maintenance of open-source software used for fundamental research in mathematics, statistics and theoretical physics. We hope that this group will consider what support is currently available, and whether there are projects that the Simons Foundation could undertake that would add significantly to the usefulness of computational tools for basic research. Modes of support that duplicate or marginally improve on support that is already available through the universities or the federal government will not be of interest to the foundation. Questions of software that is primarily educational in nature may be useful as a comparison, but are not of primary interest.  The scale of foundation support will depend upon what is needed and on the potential scientific benefit, but could be substantial, perhaps up to several million dollars per year.
Current modes of funding for research software in mathematics, statistics and physics differ very significantly. There may be correspondingly great differences in what the foundation might accomplish in these areas. We hope that the round table members will be able to help the foundation understand the current landscape  (what are the needs, what is available, whether it is useful, how it is supported) both in general and across the different disciplines, and will help us think creatively about new possibilities.
I flew across country to this the meeting, where we spent the day discussing ways in which "several million dollars per year" could revolutionize "the development, deployment and maintenance of open-source software used for fundamental research in mathematics...".

In the afternoon Jim Simons arrived, and shook our hands. He then lectured us with some anecdotes, didn't listen to what we had to say, and didn't seem to understand open source software. I was frustrated watching how he treated the other participants, so I didn't say a word to him. I feel bad for failing to express myself.

The Decision

In the backroom during a coffee break, David Eisenbud told me that it had already been decided that they were going to just fund Magma by making it freely available to all academics in North America. WTF? I explained to David that Magma is closed source and that not only does funding Magma not help open source software like Sage, it actively hurts it. A huge motivation for people to contribute to Sage is that they do not have access to Magma (which was very expensive).

I wandered out of that meeting in a daze; things had gone so differently than I had expected. How could a goal to "facilitate the development, deployment and maintenance of open-source software... perhaps up to several million dollars per year" result in a decision that would make things possibly much worse for open source software?

That day I started thinking about creating what would become SageMathCloud. The engineering work needed to make Sage accessible to a wider audience wasn't going to happen without substantial funding (I had put years of my life into this problem but it's really hard, and I couldn't do it by myself). At least I could try to make it so people don't have to install Sage (which is very difficult). I also hoped a commercial entity could provide a more sustainable source of funding for open source mathematics software. Three years later, the net result of me starting SageMathCloud and spending almost every waking moment on it is that I've gone from having many grants to not, and SageMathCloud itself is losing money. But I remain cautiously optimistic and forge on...

We will not fund Sage

Prompted by numerous messages recently from people, I wrote to David Eisenbud this week. He suggested I write to Yuri Schinkel, who is the current director of the Simons Foundation:
Dear William,
Before I joined the foundation, there was a meeting conducted by David Eisenbud to discuss possible projects in this area, including Sage.
After that meeting it was decided that the foundation would support Magma.
Please keep me in the loop regarding developments at Sage, but I regret that we will not fund Sage at this time.
Best regards, Yuri
The Simons Foundation, the NSF, or any other foundation does not owe the Sage project anything. Sage is used by a lot of people for free, who together have their research and teaching supported by hundreds of millions of dollars in NSF grants. Meanwhile the Sage project barely hobbles along. I meet people who have fantastic development or documentations projects for Sage that they can't do because they are far too busy with their fulltime teaching jobs. More funding would have a massive impact. It's only fair that the US mathematical community is at least aware of a missed opportunity.
Funding in Europe for open source math software is much better.

Hacker News discussion

by William Stein (noreply@blogger.com) at September 05, 2015 04:47 PM

Vince Knight

Picking a good Vainglory jungler with game theory and sagemath

I’ve recently been playing a really cool video game: Vainglory. This is described as a MOBA which I must admit I had never heard off until this year when my students mentioned it to me, but basically it’s an online multi player game in which players form two teams of 6 heroes and fight each other. The choice of the heroes is very important as the composition of a team can make or break a match. This seems to have a bit of a cult following (so no doubt just like for my post about clash of clans I might annoy people again) and there is a great wiki that gives guides for the play of each player. In this post I’ll describe using Python to scrape that wiki to get data that feeds in to a game theoretic model which I then analyse using Sagemath to give some insight about the choice of hero.

Here’s the map where this all takes place:


So first of all, my understanding is that there are generally three types of playing strategy in the game:

  • Lane: a hero that occupies and tries to take over the main route between the two bases.
  • Jungle: a hero that goes ‘off road’ and kills monsters, gets gold etc…
  • Roam: a hero who roams in between the two and whose main job is to support the other two players.

My personal strategy is to pick a roamer/protector: Ardan (pic below),


I generally help out the jungler in my team and try my best to not be a liability by dying.

The wiki has a bunch of information for players. If you google something like ‘vainglory best strategy’ it comes up. If you look up each hero you get a collection of guides ranked by votes each with all sorts of information which includes the where each and every other hero sits on a threat level (from 1 to 10). Here is the threat meter for Ardan from the top guide:

Threat for Ardan

So from that guide it looks like if your opponent is going to be isolated with Ardan then you should pick HERO. In some guides the threat meter does not list all the heros. This is particularly important as it’s these threat meters that I’ve used as a source of data for how good a given hero is against other heros.

This is where the keener player/reader will note that the threat meter only describes the threat to a single player and not any information about how this fits within a team dynamic. This is an important admission on my part: as indicated by the title of this post aims to use data and game theory to give an indication as to how to choose heros for isolated combat against other single heros. So one application of this is choosing a jungler/laner when you expect to go up against another team that is playing a single jungler/laner.

Scraping the data

First things first: I used Python with the BeautifulSoup and requests library. For example here is how I got the lists of all the heroes (and the url to their own respective page on the wiki):

>>> page = requests.get('http://www.vaingloryfire.com/vainglory/wiki/heroes')
>>> soup = BeautifulSoup(page.text, 'html.parser')
>>> root = '/vainglory/wiki/heroes'
>>> urls = [link.get('href') for link in soup.find_all('a')]
>>> heroes = {url[len(root) + 1:]:url for url in urls[2:] if url.startswith(root + '/')}
>>> del heroes['skye'] # Removing skye as she is brand new
{u'adagio': u'/vainglory/wiki/heroes/adagio',
 u'ardan': u'/vainglory/wiki/heroes/ardan',
 u'catherine': u'/vainglory/wiki/heroes/catherine',
 u'celeste': u'/vainglory/wiki/heroes/celeste',
 u'fortress': u'/vainglory/wiki/heroes/fortress',
 u'glaive': u'/vainglory/wiki/heroes/glaive',
 u'joule': u'/vainglory/wiki/heroes/joule',
 u'koshka': u'/vainglory/wiki/heroes/koshka',
 u'krul': u'/vainglory/wiki/heroes/krul',
 u'petal': u'/vainglory/wiki/heroes/petal',
 u'ringo': u'/vainglory/wiki/heroes/ringo',
 u'rona': u'/vainglory/wiki/heroes/rona',
 u'saw': u'/vainglory/wiki/heroes/saw',
 u'skaarf': u'/vainglory/wiki/heroes/skaarf',
 u'taka': u'/vainglory/wiki/heroes/taka',
 u'vox': u'/vainglory/wiki/heroes/vox'}

(Note there that I’m removing a brand new hero: Skye as she was released pretty much at the same time as I was writing this post.)

You can see the JuPyTer notebook which shows the code. The main technicality is that I only scraped guides from the front page for each hero. As I’ll describe later, I ran my analysis taking the average threats for a variety of cases: only taking the first guide, only taking the first 2 guides, the first 3 guides etc…

Here for example is the threats data for Adagio if you only look at this first guide:

[0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 4, 0, 0, 7, 0, 0]

Cross referencing that with the order given by the list of heroes above we see that Skaarf ranks a 7 on the threat meter to Adagio, and Ringo and Joule a 4. All the 0s are what I’ve decided to do when a threat meter does not include a given hero: indicating that that hero is not a threat to that hero. I don’t really like this as a solution but it’s probably the least worst way to deal with it (if anyone has a better way of handling this please let me know in the comments).

Here is the threats data for Krul:

[6, 3, 4, 3, 6, 4, 3, 7, 5, 5, 4, 0, 6, 6, 5, 0]

We see that in this case the only heroes that pose no threat to Krul are Fortress and Rona. Thus if your opponent is playing those heroes Krul is a best response.

As will be described in the next section, we need to build up a matrix of these rows which basically shows how well a given hero does against others. Here is the matrix of this when considering the row players and taking the opposite of the threats when using just the top guide:

If you consider a column (that corresponds to a hero) of that matrix, the row player aims to find the row that gives the highest score, which because we’ve taken the opposite of the threat score corresponds to minimising the threat posed by the column hero. This is in essence a risk averse approach, at the very end I’ll comment on what happens to the results when players aim to maximise the threat they pose.

Now that I’ve described the data (you can find all the data written to specific csv files here) I’ll go on to talk about the game theory used to try and see what the equilibrium choice of strategies should/could be.

Game theoretic analysis

All of this has been done using Sagemath, a great open source mathematics package that offers an alternative to Maple, Mathematica etc…

If you’re not familiar with game theory, this video might help (it shows the basics of game theory and how Sagemath can be used to find Nash equilibria):

Before talking about equilibria let’s just look at best response dynamics.

Using Sage we can first of all build up the normal form game for a given number of guides:

sage: def build_game(row_player_file, col_player_file):
....:    """Import the bi matrices and create the game object"""
....:    bi_matrices = []
....:    for fle in [row_player_file, col_player_file]:
....:        f = open(fle, 'r')
....:        csvrdr = csv.reader(f)
....:        bi_matrices.append(-matrix([[float(ele) for ele in row] for row in csvrdr]))
....:        f.close()
....:    return NormalFormGame(bi_matrices)
sage: g = build_game("A-01.csv", "B-01.csv")

Using this and the best_response method on Sagemath NormalFormGames we can build up all the best responses (according to a given number of guides) go each player. The cool thing is that Sagemath has some awesome graph theory written in there so we can transform that in to a nice picture (again: all the code for this can be found here):

best response graph for 1st guide

That plot confirms what we have seen earlier, we see that Krul is a best response to Fortress or Rona. Sadly, because there are so many zeros when just using the first guide, there are a bunch of heros that are not considered a threat to any of the players so they have multiple best responses and our graph is messy.

Here is the best response graph when taking the mean threats over all front page guides:

best response graph for all guides

Note that Game Theory assumes that everyone know that everyone know that everyone knows… all this. So for example if two players both player Adagio we are at an equilibrium. However if one player plays Saw then the graph indicates that the opponent should play Koshka, which means that the first player should then deviate and play Fortress which is then also an equilibrium (bot players are playing best responses to each other).

From here on I will continue the analysis using the average utility from all the guides (I’ll come back to this at the end).

So we can use Sagemath to compute all the equilibria for us. A Nash equilibria need not be a pure strategy and so will at times be a probability vector indicating how players should randomly pick a hero. Here for example is the 4th equilibrium computed by Sagemath:

sage: g.obtain_nash(algorithm='lrs')[3]
[(0, 0, 0, 0, 0, 0, 0, 0, 3947/17781, 0, 3194/17781, 0, 8795/17781, 0, 0, 615/5927),
 (0, 0, 0, 0, 0, 0, 0, 0, 3947/17781, 0, 3194/17781, 0, 8795/17781, 0, 0, 615/5927)]

This particular equilibria has both players playing a mix of: Fortress, Glaive, Petal and Koshka.

Here is the mean probability distribution for both players, while the particular values should be ignored what is of interest is the heroes that are not played at all. In essence these heroes, accross all the equilibria are not deemed playable:

ne graph for all

We see that this confirms how the previous graph was colored showing the heroes that should be played in blue.

Note that the number of guides and the reliability of all this has a huge effect of the conclusions made. Here are two gifs that show the effect of the number of guides used:

best response dynamics animation

ne graph animation

and here is a plot of the number of equilibria for each guide:

number of equilibria

Up until now all the results are for when players aim to minimise the threat posed to them. In the next section I’ll invert that (python wise it’s a minor swapping around of some inputs) and consider the situation where you want to pick a hero that is aims to be the most threatening.

Seeking to be a threat

First of all here is the best response graph:

best response graph for all

Here is the average of the NE:

best response graph for all

Those 3 players have certainly been able to rip through me on more than one occasion…

Finally here are the Nash equilibria for when a threatening player (plotted in black) is playing against a threat averse player (plotted in grey):

best response graph for all


The main thing that needs to be highlighted before concluding is that this analysis has two weaknesses:

  • The data: what comes out of mathematical models is only as good as what goes in. Scraping the wiki data is a cool thing to do (from a Python point of view) but I’m blindly grabbing guides that might have poor information/opinions in them. This is worth remembering. If someone where to come up with their own threat/performance measures then this work could just be used on that. Ultimately the data available here is better than no data.
  • I am not taking in to account team dynamics. I’m just looking at perceived threats from one hero to another. There are mathematical approaches that could be used to find the best combination of teams and I might get to that in other post one day. Nonetheless this has been a fun application of game theory and still has value I believe.

So to conclude, basing things on the data available to me, I’d suggest that (when both players are acting in a risk averse way) the choice of heros for an isolated job like jungling and/or laneing is in fact reduced to a set from:

If you and your opponent aim to be threatening, the choice is from:

Finally if you aim to be threatening, playing against a player aiming to be risk averse the choice is from:

and vice versa:

(Interestingly for this last type of game there were in general just 1 equilibrium.)

Based on all of this, I would suggest (looking across all of that summary) that (disclaimer: based on the wiki threat data) the best Vainglory hero is Glaive. Again though, this does not take in to account any of the very important team dynamics. I plan to keep on being a protector with Ardan and just doing my best to stay alive…

Another point is that this shows that vainglory is perhaps not immediately balanced. A perfectly balanced game (like Rock Paper Scissor for example) has a Nash Equilibria that evenly plays all strategies:

sage: g = game_theory.normal_form_games.RPS()
sage: g.obtain_nash()
[[(1/3, 1/3, 1/3), (1/3, 1/3, 1/3)]]

Please do take a look at all the code/data at this repository.

This was a fun application of mathematical modelling, I also learnt how to scrape with BeautifulSoup but I mainly look forward to using this in my game theory class this year. I might even suggest we spend 25 minutes of one class having a game on the big screen assuming there are 5 players of Vainglory in my class.

September 05, 2015 12:00 AM

September 01, 2015

William Stein

React, Flux, RethinkDB and SageMathCloud -- Summer 2015 update

I've been using databases and doing web development for over 20 years, and I've never really loved any database before and definitely didn't love any web development frameworks either. That all changed for me this summer...


SageMathCloud is a web application in which you collaboratively use Python, LaTeX, Markdown, Sage worksheets (sophisticated mathematics), task lists, R, Jupyter Notebooks, manage courses, write C programs, make chatrooms, and more. It is hosted on Google Compute Engine, but is also entirely open source and there is a pre-made Virtual Machine that you can download. A project in SMC is a Linux account, with resources constrained using cgroups and quotas. Many SMC users can collaborate on the same project, and have equal privileges in that project. Interaction with all file types (including Jupyter notebooks, task lists and course managements) is synchronized in realtime, like Google docs. There is also a global notifications feed that shows all editing activity on all files in all projects on which the user collaborates, which is a sort of highly technical version of Facebook's feed.

Rewrite motivation

I originally wrote the SageMathCloud frontend using progressive-refinement jQuery (no third-party framework beyond that) and the Cassandra database. These were reasonable choices when I started. There are much better approaches now, which are critical to dramatically improving the user experience with SMC, and also growing the developer base. So far SMC has had no nontrivial outside contributions, probably due to the difficulty of understanding the code. In fact, I think nobody besides me has ever even installed SMC, despite these install notes.

We (me, Jon Lee, Nicholas Ruhland) are currently completely rewriting the entire frontend of SMC using React.js, Flux, and RethinkDB. We started this rewrite in June 2015, with Jon being supported by Google Summer of Code (2015), Nich being supported some by NSF grants from Randy Leveque and Rekha Thomas, and with me being unemployed.

Terrible funding situation

I'm living on credit cards -- I have no NSF grant support anymore, and SageMathCloud is still losing a lot of money every month, and I'm unhappy about this situation. It was either completely quit working on SMC and instead teach or consult a lot, or lose tens of thousands of dollars. I am doing the latter right now. I was very caught off guard, since this is my first summer ever to not have NSF support since I got my Ph.D. in 2000, and I didn't expect to have my grant proposals all denied (which happened in June). There is some modest Angel investment in SageMath, Inc., but I can't bring myself to burn through that money on salary, since it would run out quickly, and I don't want to have to shut down the site due to not being able to pay the hosting bill. I've failed to get any significant free hosting, due to already getting free hosting in the past, and SageMath, Inc. not being in any incubators. For example, we tried very hard to get hosting from Google, but they flatly refused for these two reasons (they gave $60K in hosting to UW/Sage project in 2012). I'm clearly having trouble transitioning from an academic to an industry funding model. But if there are enough paying customers by January 2016, things will turn around.

Jon, Nich, and I have been working on this rewrite for three months, and hope to finish it by the end of September, when Jon and Nich will become busy with classes again. However, it seems unlikely we'll be able to finish at the current rate. Fortunately, I don't start teaching fulltime again until January, and we put a lot of work into doing a release in mid-August that fully uses RethinkDB and partly uses React.js, so that we can finish the second stage of the rewrite iteratively, without any major technical surprises.


Cassandra is an excellent database for many applications, but it is not the right database for SMC and I'm making no further use of Cassandra. SMC is a realtime application that does a lot more reading than writing to the database, and SMC greatly benefits from realtime push updates from the database. I've tried quite hard in the past to build an appropriate architecture for SMC on top of Cassandra, but it is the wrong tool for the job. RethinkDB scales up linearly (with sharding and replication), and has high availability and automatic failover as of version 2.1.2. See https://github.com/rethinkdb/rethinkdb/issues/4678 for my painful path to ensuring RethinkDB actually works for me (the RethinkDB developers are incredibly helpful!).


I learned about React.js first from some "random podcast", then got more interested in it when Chris Swenson gave a demo at a Sage Days workshop in San Diego in May 2015. React (+Flux) is a web development framework that actually has solid ideas behind it, backed by an implementation that has been optimized and tested by a highly nontrivial real world application: namely the Facebook website. Even if I were to have the idea of React, implementing in a way that is actually usable would be difficult. The key idea of React.js is that -- surprisingly -- it is possible to write efficient client-side code that describes how to render the application purely as a function of its state.

React is different than jQuery. With jQuery, you write lots of code explaining how to transform the user interface of your application from one complicated state (that you might never have anticipated happening) to another complicated state. When using React.js you don't write code about how your application's visible state changes -- instead you write code to answer the question: "given this state, what should the application look like". For me, it's a game changer. This is like what one does when writing video games; the innovation is that some people at Facebook figured out how to practically program this way in a client side web browser application, then tuned their implementation based on huge amounts of real world data (Facebook has users). Oh, and they open sourced the result and ran several conferences explaining React.

React.js reminds me of when Andrew Wiles proved Fermat's Last Theorem in the mid 1990s. Wiles (and Ken Ribet) had genuine new ideas, which dramatically reshaped the landscape of number theory. The best number theorists quickly realized this and adopted to the new world, pushing the envelope of Wiles work far beyond what I expected could happen. Other people pretended like Wiles didn't exist and continued studying Fibonnaci numbers. I browsed the web development section of Barnes and Noble last night and there were dozens of books on jQuery and zero on React.js. I feel for anybody who tries to learn client-side web development by reading books at Barnes and Noble.

IPython/Jupyter and PhosphorJS

I recently met with Fernando Perez, who founded IPython/Jupyter. He seemed to tell me that currently 9 people are working fulltime on rewriting the Jupyter web notebook using the PhosphorJS framework. I tried to understand PhosphorJS based on the github page, but couldn't, except to deduce that it is mostly the work of one person from Bloomberg/Continuum Analytics. Fernando told me that they chose PhosphorJS since it very fast, and that their main motivation is to (1) make Jupyter better use their huge high-resolution monitors on their new institute at Berkeley, and (2) make it easier for developers like me to integrate/extend Jupyter into their applications. I don't understand (2), because PhosphorJS is perhaps the least popular web framework I've ever heard of (is it a web framework -- I can't tell?). I pushed Fernando to explain why they made that design choice, but didn't really understand the answer, except that they had spent a lot of time investigating alternatives (like React first). I'm intimidated by their resources and concerned that I'm making the wrong choice; however, I just can't understand why they have made what seems to me to be the wrong choice. I hope to understand more at the joint Sage/Jupyter Days 70 that we are organizing together in Berkeley, CA in November. (Edit: see https://github.com/ipython/ipython/issues/8239 for a discussion of why IPython/Jupyter uses PhosphorJS.)

Tables and RethinkDB

Our rewrite of SMC is built on Tables, Flux and React. Tables are client-side technology I wrote inspired by Facebook's GraphQL/Relay technology (and Meteor, Firebase, etc.); they synchronize data between clients and the backend database in realtime. Tables are defined by a JSON schema file, which specifies the fields in the table, and explains what get and set queries are allowed. A table is a subset of a much larger table in the database, with the subset defined by conditions that are relative to the user making the query. For example, the projects table has one entry for each project that the user is a collaborator on.

Tables are automatically synchronized between the user and the database whenever the database changes, using RethinkDB changefeeds. RethinkDB's innovation is to build realtime updates -- triggered when the result of a query to the database changes -- directly into the database at the lowest level. Of course it is possible to build something that looks the same from the outside using either a message queue (say using RabbitMQ or ZeroMQ), or by watching the replication stream from the database and triggering actions based on that (like Meteor does using MongoDB). RethinkDB's approach seems better to me, putting the abstraction at the right level. That said, based on mailing list traffic, searches, etc., it seems that very, very few people get RethinkDB yet. Also, despite years of development, RethinkDB only became "production ready" a few months ago, and only got automatic failover a few weeks ago. That said, after ironing out some kinks, I'm now using it with heavy traffic in production and it works very well.


Once data is automatically synchronized between the database and web browsers in realtime, we can build everything else on top of this. Facebook also introduced an architecture pattern that they call Flux, which works well with React. It's very different than MVC-style two-way binding frameworks, where objects are directly linked to UI elements, with an object changing causing the UI element to change and vice versa. In SMC each major part of the system has two objects associated to it: Actions and Stores. We think of them in terms of the classical CQRS pattern -- command query responsibility segregation. Actions are commands -- they are Javascript "functions" that get stuff done, but they do not return values; instead, they impact the state of the store. The store has functions that allow one to query for the state of the store, but they do not change the state of the store. The store functions must only be functions of the internal state of the store and nothing else. They might cache their results and format their output to be very convenient for rendering. But that's it.

Actions usually cause the corresponding store (or stores) to change. When a store changes, it emit a change event, which causes any React components that depend on the store to be updated, which in many cases means they are re-rendered. There are optimizations one can introduce to reduce the amount of re-rendering, which if one isn't careful leads to subtle bugs; pretty much the only subtle React UI bugs one hits are caused by such optimizations. When the UI re-renders, the user sees their view of the world change. The user then clicks buttons, types, etc., which triggers actions, which in turn update stores (and tables, hence propogating changes to the ultimate source of truth, which is the RethinkDB database). As stores update, the UI again updates, etc.


So far, we have completely (re-)written the project listing, file manager, help/status page, new file page, project log, file finder, project settings, course management system, account settings, billing, project upgrade system, and file use notifications using React, Flux, and Tables, and the result works well. Bugs are much easier to fix, and it is easy (possible?) to understand the state of the system, since it is defined by the state of the database and the corresponding client-side stores. We've completely rethought everything about the UI in doing the rewrite of the above components, and it has taken several months. Also, as mentioned above, I completely rewrote most of the backend to use RethinkDB instead of Cassandra. There were also the weeks of misery for me after we made the switch over. Even after weeks of thinking/testing/wondering "what could go wrong?", we found out all kinds of surprising little things within hours of pushing everything into production, which took more than a week of sleep deprived days to sort out.

What's left? We have to rewrite the file editor tabs system, the project tabs system, and all the applications (except course management): editing text files using Codemirror, task lists (which are suprisingly complicated!), color xterm terminals, Jupyter notebooks (which will still use an iframe for the notebook itself), Sage worksheets (with complicated html output embedded in codemirror), compressed file de-archiver, the LaTeX editor, the wiki and markdown editors, and file chat. We hope to find a clean way to abstract away the various SMC applications as plugins, so that other people can easily write their own applications/plugins that will run inside of SMC. There will be a rich collection of example plugins to build on, namely the ones listed above, which are all driven by critical-to-us real world applications.

Discussion about this blog post on Hacker News.

by William Stein (noreply@blogger.com) at September 01, 2015 07:57 AM

August 28, 2015

Vince Knight

Natural language processing of new jokes from 2015

This is a brief update to a previous post: “Python, natural language processing and predicting funny”. In that post I carried out some basic natural language processing with Python to predict whether or not a joke is funny. In this post I just update that with some more data from this year’s Edinburgh Fringe festival.

Take a look at the ipython notebook which shows graphics and outputs of all the jokes. Interestingly this year’s winning joke is not deemed funny by the basic model :) but overall was 60% right this year (which is pretty good compared to last year).

Here is a summary plot of the classifiers for different thresholds of ‘funny’:

The corresponding plot this year (with the new data):

Take a look at the notebook file and by all means grab the csv file to play (but do let me know how you get on :)).

August 28, 2015 12:00 AM

August 22, 2015

Benjamin Hackl

Google Summer of Code 2015: Conclusion

The “Google Summer of Code 2015” program has ended yesterday, on the 21. of August at 19.00 UTC. This blog entry shall provide a short wrap-up of our GSoC project.

The aim of our project was to implement a basic framework that enables us to do computations with asymptotic expressions in SageMath — and I am very happy to say that we very much succeeded to do so. An overview of all our developments can be found at meta ticket #17601.

Although we did not really follow the timeline suggested in my original proposal (mainly because the implementation of the Asymptotic Ring took way longer than originally anticipated) we managed to implement the majority of ideas from my proposal — with the most important part being that our current prototype is stable. In particular, this means that we do not expect to make major design changes at this point. Every detail of our design is well-discussed and can be explained.

Of course, our “Asymptotic Expressions” project is far from finished, and we will continue to extend the functionality of our framework. For example, although working with exponential and logarithmic terms is currently possible, it is not very convenient because the $\log$, $\exp$, and power functions are not fully implemented. Furthermore, it would be interesting to investigate the performance-gain obtained by cythonizing the central parts of this framework (e.g. parts of the MutablePoset…) — and so on…

To conclude, I want to thank Daniel Krenn for his hard work and helpful advice als my mentor, as well as the SageMath community for giving me the opportunity to work on this project within the Google Summer of Code program! :-)

by behackl at August 22, 2015 11:34 AM

August 17, 2015

Benjamin Hackl

Asymptotic Expressions: Current Developments

Since my last blog entry on the status of our implementation of Asymptotic Expressions in SageMath quite a lot of improvements have happened. Essentially, all the pieces required in order to have a basic working implementation of multivariate asymptotics are there. The remaining tasks within my Google Summer of Code project are:

  • Polish the documentation of our minimal prototype, which consists of #17716 and the respective dependencies. Afterwards, we will set this to needs_review.
  • Open a ticket for the multivariate asymptotic ring and put together everything that we have written so far there.

In this blog post I want to give some more examples of what can be done with our implementation right now and what we would like to be able to handle in the future.

Status Quo

After I wrote my last blog entry, we introduced a central idea/interface to our project: short notations. By using the short notation factory for growth groups (introduced in #18930) it becomes very simple to construct the desired growth group. Essentially, monomial growth groups (cf. #17600), i.e. groups that contain elements of the form

(for a fixed variable and powers from some base ring, e.g. the Integer Ring or even the Rational Field) are represented by
, where the base ring is also specified via its shortened name. The short notation factory then enables us to do the following:

sage: from sage.groups.asymptotic_growth_group import GrowthGroup
sage: G = GrowthGroup('x^ZZ'); G
Growth Group x^ZZ
sage: G.an_element()
sage: G = GrowthGroup('x^QQ'); G
Growth Group x^QQ
sage: G.an_element()

Naturally, this interface carries over to the generation of asymptotic rings: instead of the (slightly dubious)

keyword advertised in my last blog entry, we can now actually construct the growth group by specifying the respective growth group via its short representation:

sage: R.<x> = AsymptoticRing('x^ZZ', QQ); R
Asymptotic Ring <x^ZZ> over Rational Field
sage: (x^2 + O(x))^50
x^100 + O(x^99)

Recently, we also implemented another type of growth group: exponential growth groups (see #19028). These groups represent elements of the form

where the base is from some multiplicative group. For example, we could do the following:

sage: G = GrowthGroup('QQ^x'); G
Growth Group QQ^x
sage: G.an_element()
sage: G(2^x) * G(3^x)
sage: G(5^x) * G((1/7)^x)

Note: unfortunately, we did not implement a function that allows taking some element from some growth group (e.g.

from a monomial growth group) as the variable in an exponential growth group yet. Implementing some way to “change” between growth groups by taking the log or the exponential function is one of our next steps.

We also made this short notation a central interface for working with cartesian products. This is implemented in #18587. For example, this allows to construct growth groups containing elements like $2^x \sqrt[5]{x^2} \log(x)^2$:

sage: G = GrowthGroup('QQ^x * x^QQ * log(x)^ZZ'); G
Growth Group QQ^x * x^QQ * log(x)^ZZ
sage: G.an_element()
(1/2)^x * x^(1/2) * log(x)
sage: G(2^x * x^(2/5) * log(x)^2)
2^x * x^(2/5) * log(x)^2

Simple parsing from the symbolic ring (and from strings) is implemented. Like I have written above, operations like

are one of the next steps on our roadmap.

Further Steps

Of course, having an easy way to generate growth groups (and thus also asymptotic rings) is nice — however, it would be even better if the process of finding the correct parent would be even more automated. Unfortunately, this requires some non-trivial effort regarding the pushout construction — which will certainly not happen within the GSoC project.

As soon as we have an efficient way to “switch” between factors of a growth group (e.g. by taking the logarithm or the exponential function), this has to be carried over up to the asymptotic ring. Operations like

sage: 2^(x^2 + O(x))
2^(x^2) * 2^(O(x))

where the output could also be

2^(x^2) * O(x^g)
, where $g$ is determined by

Division of asymptotic expressions can be realized with just about the same idea, for example:

\[ \frac{1}{x^2 + O(x)} = \frac{1}{x^2} \frac{1}{1 + O(1/x)} = x^{-2} + O(x^{-3}), \]

and so on. If an infinite series occurs, it will have to be cut using an $O$-Term, most likely somehow depending on

as well.

Ultimately, we would also like to incorporate, for example, Stirling’s approximation of the factorial such that we could do something like

sage: n.factorial()
sqrt(2*pi) * e^(n*log(n)) * (1/e)^n * n^(1/2) + ...

which then can be used to obtain asymptotic expansions of binomial coefficients like $\binom{2n}{n}$:

sage: (2*n).factorial() / (n.factorial()^2)
1/sqrt(pi) * 4^n * n^(-1/2) + ...

As you can see, there is still a lot of work within our “Asymptotic Expressions” project — nevertheless, with the minimal working prototype and the ability to create cartesian products of growth groups, the fundament for all of this is already implemented! 😉

by behackl at August 17, 2015 07:15 AM

August 16, 2015

Michele Borassi

Conclusion of the Main Part of the Project

In this post, I will summarize the results obtained with the inclusion in Sage of Boost and igraph libraries. This was the main part of my Google Summer of Code project, and it was completed yesterday, when ticket 19003 was closed.

We have increased the number of graph algorithms available in Sage from 66 to 98 (according to the list used in the initial comparison of the graph libraries [1]). Furthermore, we decreased the running-time of several Sage algorithms: in some cases, we have been able to improve the asymptotic running-time, obtaining up to 10000x improvements in our tests. Finally, during the inclusion of external algorithms, we have refactored and cleaned some of Sage source code, like the shortest path routines: we have standardized the input and the output of 15 routines related to shortest paths, and we have removed duplicate code as much as possible.

More specifically, the first part of the project was the inclusion of Boost graph library: since the library is only available in C++, we had to develop an interface. This interface lets us convert easily a Sage graph into a Boost graph, and run algorithms on the converted graph. Then, we have written routines to re-translate the output into a Sage-readable format: this way, the complicated Boost library is "hidden", and users can interact with it as they do with Sage. In particular, we have interfaced the following algorithms:
  • Edge connectivity (trac.sagemath.org/ticket/18564);
  • Clustering coefficient (trac.sagemath.org/ticket/18811);
  • Cuthill-McKee and King vertex orderings (trac.sagemath.org/ticket/18876);
  • Minimum spanning tree (trac.sagemath.org/ticket/18910);
  • Dijkstra, Bellman-Ford, Johnson shortest paths (trac.sagemath.org/ticket/18931);
All these algorithms were either not available in Sage, or quite slow, compared to the Boost routines. As far as we know, Boost does not offer other algorithms that improve Sage algorithms: however, if such algorithms are developed in the future, it will be very easy to include them, using the new interface.

In the second part of the project, we included igraph: since this library already offers a Python interface, we decided to include it as an optional package (before it becomes a standard package, at least an year should pass [2]). To install the package, it is enough to type the following instruction from the Sage root folder:

sage -i igraph        # To install the igraph C core
sage -i python_igraph # To install the Python interface

Then, we can easily interact with igraph: for a list of available routines, it is enough to type "igraph." and click tab twice. To convert a Sage graph g_sage into an igraph graph it is enough to type g_igraph = g_sage.igraph_graph(), while a Sage graph can be instantiated from an igraph graph through g_sage=Graph(g_igraph) or g_sage=DiGraph(g_igraph). This way, all igraph algorithms are now available in Sage.

Furthermore, we have included the igraph maximum flow algoritm inside the Sage corresponding function, obtaining significant improvements (for more information and benchmarks, we refer to ticket 19003 [3]).

In conclusion, I think the project reached its main goal, the original plan was followed very closely, and we have been able to overcome all problems.

Before closing this post, I would like to thank many people that helped me with great advices, and who provided great solutions to all the problems I faced. First of all, my mentor David Coudert: he always answered very fast to all my queries, and he gave me great suggestions to improve the quality of the code I wrote. Then, a very big help came from Nathann Cohen, who often cooperated with David in reviewing my code and proposing new solutions. Moreover, I have to thank Martin Cross, who gave me good suggestions with Boost graph library, and Volker Braun, who closed all my ticket. Finally, I have to thank the whole Sage community for giving me this great opportunity!

[1] https://docs.google.com/spreadsheets/d/1Iu1hkQtRn9J-sgfZbQTu2RoXzyjoMEWP5-cm3nAwnWE/edit?usp=sharing
[2] http://doc.sagemath.org/html/en/developer/coding_in_other.html
[3] http://trac.sagemath.org/ticket/19003

by Michele Borassi (noreply@blogger.com) at August 16, 2015 06:59 AM

August 09, 2015

Vince Knight

Why I am a paying member of cloud.sagemath

If you are not familiar with Sagemath it is a free open source mathematics package that does simple things like expand algebraic expressions as well as far more complex things (optimisation, graph theory, combinatorics, game theory etc…). Cloud.sagemath is a truly amazing tool not just for Sage bu for scientific computation in general and it’s free. Completely 100% free. In this post I’ll explain why I pay for it.

A while ago, a colleague and I were having a chat about the fact that our site Maple license hadn’t been renewed fast enough (or something similar to that). My colleague was fairly annoyed by this saying something like:

‘We are kind of like professional athletes, if I played soccer at a professional club I would have the best facilities available to me. There would not be a question of me having the best boots.’

Now I don’t think we ever finished this conversation (or at least I don’t really remember what I said) but this is something that’s stayed with me for quite a while.

First of all:

I think there are probably a very large proportion of professional soccer players who do not play at the very top level and so do not enjoy having access to the very best facilities (I certainly wouldn’t consider myself the Ronaldo of mathematics…).


Mathematicians are (in some ways) way cooler than soccer players. We are somewhat like magicians, in the past we have not needed much more than a pencil and some paper to work our craft. Whilst a chemist/physicist/medical research needs a lab and/or other things we can pretty much work just with a whiteboard.

We are basically magicians. We can make something from nothing.

Since moving to open source software for all my research and teaching this is certainly how I’ve felt. Before discovering open source tools I needed to make sure I had the correct licence or otherwise before I could work but this is no longer the case. I just need a very basic computer (I bought a thinkpad for £60 the other day!) and I am just as powerful as I could want to be.

This is even more true with cloud.sagemath. Anyone can use a variety of scientific computing tools for no cost whatsoever (not even a cost associated with the time spent installing software): it just works. I have used this to work on sage source code with students, carry out research and also to deliver presentations: it’s awesome.

So, why do I pay $7 month to use it?

Firstly because it gives me the ability to move some projects to servers that are supposedly more robust. I have no doubt that they are more robust but in all honesty I can’t say I’ve seen problems with the ‘less’ robust servers (150 of my students used them last year and will be doing so again in the Autumn).

The main reason I pay to use cloud.sagemath is because I can afford to.

This was put in very clear terms to me during the organisation of DjangoCon Europe. The principle at Python conferences is that everyone pays to attend. This in turn ensures that funds are available for people who cannot afford to pay to attend.

I am in a lucky enough financial position that for about the price of two fancy cups of coffee a month I can help support an absolutely amazing project that helps everyone and anyone have the same powers a magician does. This helps (although my contribution is obviously a very small part of it) ensure that students and anyone else who cannot afford to help support the project, can use Sage.

August 09, 2015 12:00 AM

August 01, 2015

Vince Knight

Simulating continuous Markov chains

In a blog post I wrote in 2013, I showed how to simulate a discrete Markov chain. In this post we’ll (written with a bit of help from Geraint Palmer) show how to do the same with a continuous chain which can be used to speedily obtain steady state distributions for models of queueing processes for example.

A continuous Markov chain is defined by a transition rate matrix which shows the rates at which transitions from 1 state to an other occur. Here is an example of a continuous Markov chain:

This has transition rate matrix \(Q\) given by:

The diagonals have negative entries, which can be interpreted as a rate of no change. To obtain the steady state probabilities \(\pi\) for this chain we can solve the following matrix equation:

if we include the fact that the sum of \(\pi\) must be 1 (so that it is indeed a probability vector) we can obtain the probabilities in Sagemath using the following:

You can run this here (just click on ‘Evaluate’):

Q = matrix(QQ, [[-3, 2, 1], [1, -5, 4], [1, 8, -9]])

This returns:

(1/4, 1/2, 1/4)

Thus, if we were to randomly observe this chain:

  • 25% of the time it would be in state 1;
  • 50% of the time it would be in state 2;
  • 25% of the time it would be in state 3.

Now, the markov chain in question means that if we’re in the first state the rate at which a change happens to go to the second state is 2 and the rate at which a change happens that goes to the third state is 1.

This is analagous to waiting at a bus stop at the first city. Buses to the second city arrive randomly 2 per hour, and buses to the third city arrive randomly 1 per hour. Everyone waiting for a bus catches the first one that arrives. So at steady state the population will be spread amongst the three cities according to \(\pi\).

Consider yourself at this bus stop. As all this is Markovian we do not care what time you arrived at the bus stop (memoryless property). You expect the bus to the second city to arrive 1/2 hours from now, with randomness, and the bus to the third city to arrive 1 hour from now, with randomness.

To simulate this we can sample two random numbers from the exponential distribution and find out which bus arrives first and ‘catch that bus’:

import random
[random.expovariate(2), random.expovariate(1)]

The above returned (for this particular instance):

[0.5003491524841699, 0.6107995795458322]

So here it’s going to take .5 hours for a bus to the second city to arrive, whereas it would take .61 hours for a bus to the third. So we would catch the bust to the second city after spending 0.5 hours at the first city.

We can use this to write a function that will take a transition rate matrix, simulate the transitions and keep track of the time spent in each state:

def sample_from_rate(rate):
    import random
    if rate == 0:
        return oo
    return random.expovariate(rate)

def simulate_cmc(Q, time, warm_up):
    Q = list(Q)  # In case a matrix is input
    state_space = range(len(Q))  # Index the state space
    time_spent = {s:0 for s in state_space}  # Set up a dictionary to keep track of time
    clock = 0  # Keep track of the clock
    current_state = 0  # First state
    while clock < time:
        # Sample the transitions
        sojourn_times = [sample_from_rate(rate) for rate in Q[current_state][:current_state]]
        sojourn_times += [oo]  # An infinite sojourn to the same state
        sojourn_times += [sample_from_rate(rate) for rate in Q[current_state][current_state + 1:]]

        # Identify the next state
        next_state = min(state_space, key=lambda x: sojourn_times[x])
        sojourn = sojourn_times[next_state]
        clock += sojourn
        if clock > warm_up:  # Keep track if past warm up time
            time_spent[current_state] += sojourn
        current_state = next_state  # Transition

    pi = [time_spent[state] / sum(time_spent.values()) for state in state_space]  # Calculate probabilities
    return pi

Here are the probabilities from the same Markov chain as above:

which gave (on one particular run):

[0.25447326473556037, 0.49567517998307603, 0.24985155528136352]

This approach was used by Geraint Palmer who is doing a PhD with Paul Harper and I. He used this to verify that calculations were being carried out correctly when he was trying to fit a model. James Campbell and I are going to try to use this to get an approximation for bigger chains that cannot be solved analytically in a reasonable amount of time. In essence the simulation of the Markov chain makes sure we spend time calculating probabilities in states that are common.

August 01, 2015 12:00 AM

July 27, 2015

Michele Borassi

Including igraph Library

In this new blog post, I would like to discuss the inclusion of igraph library inside Sage.
Up to now, I have interfaced Sagemath with Boost graph library, in order to run Boost algorithms inside Sage. Now, I want to do the same with igraph, the other major C++ graph library, which stands out because it contains 62 routines, 29 of which are not available in Sage. Moreover, igraph library is very efficient, as shown in [1] and in the previous post on library comparison.

This inclusion of igraph in Sage is quite complicated, because we have to include a new external library [2] (while in the Boost case we already had the sources). We started this procedure through ticket 18929: unfortunately, after this ticket is closed, igraph will only be an optional package, and we will have to wait one year before it becomes standard. The disadvantage of optional packages is that they must be installed before being able to use them; however, the installation is quite easy: it is enough to run Sage with option -i python_igraph.

After the installation, the usage of igraph library is very simple, because igraph already provides a Python interface, that can be used in Sage. To transform the Sagemath network g_sage into an igraph network g_igraph, it is enough to type g_igraph=g_sage.igraph_graph(), while to create a Sagemath network from an igraph network it is enough to type g_sage = Graph(g_igraph) or g_sage=DiGraph(g_igraph). After this conversion, we can use all routines offered by igraph!
For instance, if we want to create a graph through the preferential attachment model, we can do it with the Sagemath routine, or with the igraph routine:

sage: G = graphs.RandomBarabasiAlbert(100, 2)
sage: G.num_verts()
sage: G = Graph(igraph.Graph.Barabasi(100, int(2)))
sage: G.num_verts()

The result is the same (apart from randomness), but the time is very different:

sage: import igraph
sage: %timeit G = Graph(igraph.Graph.Barabasi(10000000, int(2)))
1 loops, best of 3: 46.2 s per loop

sage: G = graphs.RandomBarabasiAlbert(10000000, 2)
Stopped after 3 hours.

Otherwise, we may use igraph to generate graphs with Forest-Fire algorithm, which is not available in Sagemath:

sage: G = Graph(igraph.Graph.Forest_Fire(10, 0.1))
sage: G.edges()
[(0, 1, None), (0, 2, None), (1, 7, None), (2, 3, None), (2, 4, None), (3, 5, None), (3, 8, None), (4, 6, None), (8, 9, None)]

We may also do the converse: transform a Sage network into an igraph network and apply an igraph algorithm. For instance, we can use label propagation to find communities (a task which is not implemented in Sage):

sage: G = graphs.CompleteGraph(5)+graphs.CompleteGraph(5)
sage: G.add_edge(0,5)
sage: com = G.igraph_graph().community_label_propagation()
sage: len(com)
sage: com[0]
[0, 1, 2, 3, 4]
sage: com[1]
[5, 6, 7, 8, 9]

The algorithm found the two initial cliques as communities.

I hope that these examples are enough to show the excellent possibilities offered by igraph library, and that these features will soon be available in Sagemath!

[1] https://sites.google.com/a/imtlucca.it/borassi/unpublished-works/google-summer-of-code/library-comparison
[2] http://doc.sagemath.org/html/en/developer/packaging.html

by Michele Borassi (noreply@blogger.com) at July 27, 2015 09:07 AM

Vince Knight

A talk on computational game theory in Sagemath

Today, Cardiff University, School of Mathematics students: James Campbell, Hannah Lorrimore as well as Google Summer of Code student Tobenna P. Igwe (PhD student at the University of Liverpool) as well as I presented the current game theoretic capabilities of Sagemath.

This talk happened as part of a two day visit to see Dima Pasechnik to work on the stuff we’ve been doing and the visit was kindly supported by CoDiMa (an EPSRC funded project to support the development of GAP and Sagemath)

Here is the video of the talk:

Here is a link to the sage worksheet we used for the talk.

Here are some photos I took during the talk:

and here are some I took of us working on code afterwards:

Here is the abstract of the talk:

Game Theory is the study of rational interaction and is getting increasingly important in CS. Ability to quickly compute a solution concept for a nontrivial (non-)cooperative game helps a lot in practical and theoretic work, as well as in teaching. This talk will describe and demonstrate the game theoretic capabilities of Sagemath (http://www.sagemath.org/), a Python library, described as having the following mission: ‘Creating a viable free opensource alternative to Magma, Maple, Mathematica and Matlab’.

The talk will describe algorithms and classes that are implemented for the computation of Nash equilibria in bimatrix games. These include:

  • A support enumeration algorithm;
  • A reverse search algorithm through the lrs library;
  • The Lemke-Howson algorithm using the Gambit library (https://github.com/gambitproject/gambit).

In addition to this, demonstrations of further capabilities that are actively being developed will also be given:

  • Tests for degeneracy in games;
  • A class for extensive form games which include the use of the graph theoretic capabilities of Sage.

The following two developments which are being carried out as part of a Google Summer of Code project will also be demonstrated:

  • An implementation of the Lemke-Howson algorithm;
  • Extensions to N player games;

Demonstrations will use the (free) online tool cloud.sagemath which allows anyone with connectivity to use Sage (and solve game theoretic problems!). Cloud.sagemath also serves as a great teaching and research tool with access to not only Sage but Jupyter (Ipython) notebooks, R, LaTeX and a variety of other software tools.

The talk will concentrate on strategic non-cooperative games but matching games and characteristic function games will also be briefly discussed.

July 27, 2015 12:00 AM

July 23, 2015

Vince Knight

Using the two thirds of the average game in class

This past week I have been delighted to have a short pedagogic paper accepted for publication in MSOR Connections. The paper is entitled: “Playing Games: A Case Study in Active Learning Applied to Game Theory”. The journal is open access and you can see a pre print here. As well as describing some literature on active learning I also present some data I’ve been collecting (with the help of others) as to how people play two subsequent plays of the two thirds of the average game (and talk about another game also).

In this post I’ll briefly put up the results here as well as mention a Python library I’m working on.

If you’re not familiar with it, the two thirds of the average game asks players to guess a number between 0 and 100. The closest number to 2/3rds of the average number guessed is declared the winner.

I use this all the time in class and during outreach events. I start by asking participants to play without more explanation than the basic rules of the game. Following this, as a group we go over some simple best response dynamics that indicate that the equilibrium play for the game is for everyone to guess 0. After this explanation, everyone plays again.

Below you can see how this game has gone as a collection of all the data I’ve put together:

You will note that some participants actually increase their second guess but in general we see a possible indication (based on two data points, so obviously this is not meant to be a conclusive statement) of convergence towards the theoretic equilibria.

Here is a plot showing the relationship between the first and second guess (when removing the guesses that increase, although as you can see in the paper this does not make much difference):

The significant linear relationship between the guesses is given by:

So a good indication of what someone will guess in the second round is that it would be a third of their first round guess.

Here is some Sage code that produces the cobweb diagram assuming the following sequence represents each guess (using code by Marshall Hampton):

that plot shows the iterations of the hypothetical guesses if we were to play more rounds :)

The other thing I wanted to point at in this blog post is this twothirds library which will potentially allow anyone to analyse these games quickly. I’m still working on it but if it’s of interest please do jump in :) I have put up a Jupyter notebook demoing what it can do so far (which is almost everything but with some rough edges). If you want to try it out, download that notebook and run:

$ pip install twothirds

I hope that once the library is set up anyone who uses it could simply send over data of game plays via PR which would help update the above plots and conclusions :)

July 23, 2015 12:00 AM

July 16, 2015

Benjamin Hackl

Computing with Asymptotic Expressions

It has been quite some time since my last update on the progress of my Google Summer of Code project, which has two reasons. On the one hand, I have been busy because of the end of the semester, as well as because of the finalization of my Master’s thesis — and on the other hand, it is not very interesting to write a post on discussing and implementing rather technical details. Nevertheless, Daniel Krenn and myself have been quite busy in order to bring asymptotic expressions to SageMath. Fortunately, these efforts are starting to become quite fruitful.

In this post I want to discuss our current implementation roadmap (i.e. not only for the remaining Summer of Code, but also for the time afterwards), and give some examples for what we are currently able to do.

Strutcture and Roadmap

An overview of the entire roadmap can be found at here (trac #17601). Recall that the overall goal of this project is to bring asymptotic expressions like $2^n + n^2 \log n + O(n)$ to Sage. Our implementation (which aims to be as general and expandable as possible) tackles this problem with a three-layer approach:

  • GrowthGroups and GrowthElements (trac #17600). These elements and parents manage the growth (and just the growth!) of a summand in an asymptotic expression like above. The simplest cases are monomial and logarithmic growth groups. For example, their elements are given by $n^r$ and $\log(n)^r$ where the exponent $r$ is from some ordered ring like $\mathbb{Z}$ or $\mathbb{Q}$. Both cases (monomial and logarithmic growth groups) can be handled in the current implementation — however, growth elements like $n^2 \log n$ are intended to live in the cartesian product of a monomial and a logarithmic growth group (in the same variable). Parts of this infrastructure are already prepared (see trac #18587).
  • AsymptoticTerms and TermMonoids (trac #17715). While GrowthElements only represent the growth, AsymptoticTerms have more information: basically, they represent a summand in an asymptotic expression. There are different classes for each type of asymptotic term (e.g. ExactTerm and OTerm — with more to come). Additionally to a growth element, some types of asymptotic terms (like exact terms) also possess a coefficient.
  • AsymptoticExpression and AsymptoticRing (trac #17716). This is what we are currently working on, and we do have a running prototype! :-) The version that can be found on trac is only missing some doctests and a bit of documentation. Asymptotic expressions are the central objects within this project, and essentially they are sums of several asymptotic terms. In the background, we use a special data structure (“mutable posets“, trac #17693) in order to model the (partial) order induced by the various growth elements belonging to an asymptotic expression. This allows to perform critical operations like absorption (when an \(O\)-term absorbs “weaker” terms) efficiently and in a simple way.

The resulting minimal prototype can, in some sense, be compared to Sage’s PowerSeriesRing: however, we also allow non-integer exponents, and extending this prototype to work with multivariate expressions should not be too hard now, as the necessary infrastructure is there.

Following the finalization of the minimal prototype, there are several improvements to be made. Here are some examples:

  • Besides addition and multiplication, we also want to divide asymptotic expressions, and higher-order operations like exponentiation and taking the logarithm would be interesting as well.
  • Also, conversion from, for example, the symbolic ring is important when it comes to usability of our tools. We will implement and enhance this conversion gradually.


An asymptotic ring (over a monomial growth group with coefficients and exponents from the rational field) can be created with

sage: R.<x> = AsymptoticRing('monomial', QQ); R
Asymptotic Ring over Monomial Growth Group in x over Rational Field with coefficients from Rational Field

Note that we marked the code as experimental, meaning that you will see some warnings regarding the stability of the code. Now, as we have an asymptotic ring, we can do some calculations. For example, take $ (2\sqrt{x} + O(1))^{15}$:

sage: (2*x^(1/2) + O(x^0))^15
O(x^7) + 32768*x^(15/2)

We can also have a look at the underlying structure:

sage: expr = (x^(3/7) + 2*x^(1/5)) * (x + O(x^0)); expr
O(x^(3/7)) + 2*x^(6/5) + 1*x^(10/7)
sage: expr.poset
poset(O(x^(3/7)), 2*x^(6/5), 1*x^(10/7))
sage: print expr.poset.full_repr()
poset(O(x^(3/7)), 2*x^(6/5), 1*x^(10/7))
+-- null
|   +-- no predecessors
|   +-- successors:   O(x^(3/7))
+-- O(x^(3/7))
|   +-- predecessors:   null
|   +-- successors:   2*x^(6/5)
+-- 2*x^(6/5)
|   +-- predecessors:   O(x^(3/7))
|   +-- successors:   1*x^(10/7)
+-- 1*x^(10/7)
|   +-- predecessors:   2*x^(6/5)
|   +-- successors:   oo
+-- oo
|   +-- predecessors:   1*x^(10/7)
|   +-- no successors

As you might have noticed, the “O”-constructor that is used for the PowerSeriesRing and related structures, can also be used here. In particular, $O(\mathit{expr})$ acts exactly as expected:

sage: expr
O(x^(3/7)) + 2*x^(6/5) + 1*x^(10/7)
sage: O(expr)

Of course, the usual rules for computing with asymptotic expressions hold:

sage: O(x) + O(x)
sage: O(x) - O(x)

So far, so good. Our next step is making the multivariate growth groups usable for the AsymptoticRing and then improving the overall user interface of the ring.


by behackl at July 16, 2015 03:42 PM

July 09, 2015

Michele Borassi

New Boost Algorithms

My Google Summer of Code project is continuing, and I am currently trying to include more Boost algorithms in Sage. In this post, I will make a list of the main algorithms I'm working on.

Clustering Coefficient

If two different people have a friend in common, there is a high chance that they will become friends: this is the property that the clustering coefficient tries to capture. For instance, if I pick two random people, very probably they will not know each other, but if I pick two of my acquaintances, very probably they will know each other. In this setting, the clustering coefficient of a person is the probability that two random acquaintances of this person know each other. In order to quantify this phenomenon, we can formalize everything in terms of graphs: people are nodes and two people are connected if they are acquaintances. Hence, we define the clustering coefficient of a vertex \(v\) in a graph \(G=(V,E)\) as:
$$\frac{2|\{(x,y) \in E:x,y \in N_v\}|}{\deg(v)(\deg(v)-1)}$$ where \(N_v\) is the set of neighbors of \(v\) and \(\deg(v)\) is the number of neighbors of \(v\). This is exactly the probability that two random neighbors of \(v\) are linked with an edge.
My work has included in Sagemath the Boost algorithm to compute the clustering coefficient, which is more efficient that the previous algorithm, which was based on NetworkX:

sage: g = graphs.RandomGNM(20000,100000)
sage: %timeit g.clustering_coeff(implementation='boost')
10 loops, best of 3: 258 ms per loop
sage: %timeit g.clustering_coeff(implementation='networkx')
1 loops, best of 3: 3.99 s per loop

But Nathann did better: he implemented a clustering coefficient algorithm from scratch, using Cython, and he managed to outperform the Boost algorithm, at least when the graph is dense. Congratulations, Nathann! However, when the graph is sparse, Boost algorithm still seems to be faster.

Dominator tree

Let us consider a road network, that is, a graph where vertices are street intersections, and edges are streets. The question is: if I close an intersection, where am I still able to go, assuming I am at home?
The answer to this question can be summarized in a dominator tree. Assume that, in order to go from my home to my workplace, I can choose many different paths, but all these paths pass through the café, then they pass through the square (that is, if either the café or the square is closed, then there is no way I can go to work). In this case, in the dominator tree, the father of my workplace is the square, the father of the square is the café, and the father of the café is my home, that is also the root of the tree. More formally, given a graph \(G\), the dominator tree of \(G\) rooted at a vertex \(v\) is defined by connecting each vertex \(x\) with the last vertex \(y \neq x\) that belongs to each path from \(v\) to \(x\) (note that this vertex always exists, because \(v\) belongs to each path from \(v\) to \(x\)).
Until now, Sagemath did not have a routine to compute the dominator tree: I have been able to include the Boost algorithm. Unfortunately, due to several suggestions and improvements in the code, the ticket is not closed, yet. Hopefully, it will be closed very soon!

Cuthill-McKee ordering / King ordering

Let us consider a graph \(G=(V,E)\): a matrix \(M\) of size \(|V|\) can be associated to this graph, where \(M_{i,j}=1\) if and only if there is an edge between vertices \(i\) and \(j\).
In some cases, this matrix can have specific properties, that can be exploited for many purposes, like speeding-up algorithms. One of this properties is bandwidth, which measures how far the matrix is from a diagonal matrix: it is defined as \(\max_{M_{i,j} \neq 0}|i-j|\). A small bandwidth might help in computing several properties of the graph, like eigenvalues and eigenvectors.
Since the bandwidth depends on the order of vertices, we can try to permute them in order to obtain a smaller value: in Sage, we have a routine that performs this task. However, this routine is very slow, and it is prohibitive even for very small graphs (in any case, finding an optimal ordering is NP-hard).
Hence, researchers have developed heuristics to compute good orderings: the most important ones are Cuthill-McKee ordering and King ordering. Boost contains both routines, but Sage does not: for this reason, I would like to insert these two functions. The code is almost ready, but part of it depends on the code of the dominator tree: as soon as the dominator tree is reviewed, I will open a ticket on these two routines!

Dijkstra/Bellman-Ford/Johnson shortest paths

Let us consider again a road network. In this case, we are building a GPS software, which has to compute the shortest path between the place where we are and the destination. The textbook algorithm that performed this task is Dijkstra algorithm, which computes the distance between the starting point and any other reachable point (of course, there are more efficient algorithms involving a preprocessing, but Dijkstra is the most simple, and its running-time is asymptotically optimal). This algorithm is already implemented in Sagemath.
Let's spice things up: what if that there are some streets with negative length? For instance, we like a street so much that we are willing to drive 100km more just to pass from that street, which is 50km long. It is like that street is -50km long!
First of all, under these assumptions, a shortest path might not exist: if there is a cycle with negative length, we may drive along that cycle all the times we want, decreasing more and more the distance to the destination. At least, we have to assume that no negative cycle exists.
Even with this assumption, Dijkstra algorithm does not work, and we have to perform Bellman-Ford algorithm, which is less efficient, but more general. Now, assume that we want something more: we are trying to compute the distance between all possible pairs of vertices. The first possibility is to run Bellman-Ford algorithm \(n\) times, where \(n\) is the number of nodes in the graph. But there is a better alternative: it is possible to perform Bellman-Ford algorithm only once, and then to modify the lengths of edges, so that all lengths are positive, and shortest paths are not changed. This way, we run Dijkstra algorithm \(n\) times on this modified graph, obtaining a better running time. This is Johnson algorithm.
Both Bellman-Ford and Johnson algorithms are implemented in Boost and not in Sagemath. As soon as I manage to create weighted Boost graphs (that is, graphs where edges have a length), I will include also these two algorithm!

by Michele Borassi (noreply@blogger.com) at July 09, 2015 12:51 PM

June 25, 2015

Michele Borassi

Edge Connectivity through Boost Graph Library

After two weeks, we have managed to interface Boost and Sagemath!

However, the interface was not as simple as it seemed. The main problem we found is the genericity of Boost: almost all Boost algorithms work with several graph implementations, which differ in the data structures used to store edges and vertices. For instance, the code that implements breadth-first search works if the adjacency list of a vertex v is a vector, a list, a set, etc. This result is accomplished by using templates [1]. Unfortunately, the only way to interface Sagemath with C++ code is Cython, which is not template-friendly, yet. In particular, Cython provides genericity through fused types [2], whose support is still experimental, and which do not offer full integration with templates [3-5].

After a thorough discussion with David, Nathann, and Martin (thank you very much!), we have found a solution: for the input, we have defined a fused type "BoostGenGraph", including all Boost graph implementations, and all functions that interface Boost and Sagemath use this fused type. This way, for each algorithm, we may choose the most suitable graph implementation. For the output, whose type might be dependent on the input type, we use C++ to transform it into a "standard" type (vector, or struct).

We like this solution because it is very clean, and it allows us to exploit Boost genericity without any copy-paste. Still, there are some drawbacks:
1) Cython fused types do not allow nested calls of generic functions;
2) Boost graphs cannot be converted to Python objects: they must be defined and deleted in the same Cython function;
3) No variable can have a generic type, apart from the arguments of generic functions.

These drawbacks will be overcome as soon as Cython makes templates and generic types interact: this way, we will be able create a much stronger interface, by writing a graph backend based on Boost, so that the user might create, convert, and modify Boost graphs directly from Python. However, for the moment, we will implement all algorithms using the current interface, which already provides genericity, and which has no drawback if the only goal is to "steal" algorithms from Boost.

As a test, we have computed the edge connectivity of a graph through Boost: the code is available in ticket 18564 [6]. Since the algorithm provided by Sagemath is not optimal (it is based on linear programming), the difference in the running time is impressive, as shown by the following tests:

sage: G = graphs.RandomGNM(100,1000)
sage: %timeit G.edge_connectivity()
100 loops, best of 3: 1.42 ms per loop
sage: %timeit G.edge_connectivity(implementation="sage")
1 loops, best of 3: 11.3 s per loop

sage: G = graphs.RandomBarabasiAlbert(300,3)
sage: %timeit G.edge_connectivity(implementation="sage")
1 loops, best of 3: 9.96 s per loop
sage: %timeit G.edge_connectivity()
100 loops, best of 3: 3.33 ms per loop

Basically, on a random Erdos-Renyi graph with 100 vertices and 1000 edges, the new algorithm is 8,000 times faster, and on a random Barabasi-Albert graph with 300 nodes and average degree 3, the new algorithm is 3,000 times faster! This way, we can compute the edge connectivity of much bigger graphs, like a random Erdos-Renyi graph with 5,000 vertices and 50,000 edges:

sage: G = graphs.RandomGNM(5,000, 50,000)
sage: %timeit G.edge_connectivity()
1 loops, best of 3: 16.2 s per loop

The results obtained with this first algorithm are very promising: in the next days, we plan to interface several other algorithms, in order to improve both the number of available routines and the speed of Sagemath graph library!

[1] https://en.wikipedia.org/wiki/Template_%28C%2B%2B%29
[2] http://docs.cython.org/src/userguide/fusedtypes.html
[3] https://groups.google.com/forum/#!topic/cython-users/qQpMo3hGQqI
[4] https://groups.google.com/forum/#!searchin/cython-users/fused/cython-users/-7cHr6Iz00Y/Z8rS03P7-_4J
[5] https://groups.google.com/forum/#!searchin/cython-users/fused$20template/cython-users/-7cHr6Iz00Y/Z8rS03P7-_4J
[6] http://trac.sagemath.org/ticket/18564

by Michele Borassi (noreply@blogger.com) at June 25, 2015 12:56 AM

Vince Knight

On testing degeneracy of bi-matrix games

We (James Campbell and Vince Knight are writing this together) have been working on implementing code in Sage to test if a game is degenerate or not. In this post we’ll prove a simple result that is used in the algorithm that we are/have implemented.

Bi-Matrix games

For a general overview of these sorts of things take a look at this post from a while ago on the subject of bi-matrix games in Sage. A bi-matrix is a matrix of tuples corresponding to payoffs for a 2 player Normal Form Game. Rows represent strategies for the first player and columns represent strategies for the second player, and each tuple of the bi-matrix corresponds to a tuple of payoffs. Here is an example:

We see that if the first player plays their third row strategy and the second player their second column strategy then the first player gets a utility of 6 and the second player a utility of 1.

This can also be written as two separate matrices. A matrix \(A\) for Player 1 and \(B\) for Player 2.

Here is how this can be constructed in Sage using the NormalFormGame class:

sage: A = matrix([[3,3],[2,5],[0,6]])
sage: B = matrix([[3,2],[2,6],[3,1]])
sage: g = NormalFormGame([A, B])
sage: g
Normal Form Game with the following utilities: {(0, 1): [3, 2], (0, 0): [3, 3],
(2, 1): [6, 1], (2, 0): [0, 3], (1, 0): [2, 2], (1, 1): [5, 6]}

Currently, within Sage, we can obtain the Nash equilibria of games:

sage: g.obtain_nash()
[[(0, 1/3, 2/3), (1/3, 2/3)], [(4/5, 1/5, 0), (2/3, 1/3)], [(1, 0, 0), (1, 0)]]

We see that this game has 3 Nash equilibria. For each, we see that the supports (the number of non zero entries) of both players’ strategies are the same size. This is, in fact, a theoretical certainty when games are non degenerate.

If we modify the game slightly:

sage: A = matrix([[3,3],[2,5],[0,6]])
sage: B = matrix([[3,3],[2,6],[3,1]])
sage: g = NormalFormGame([A, B])
sage: g.obtain_nash()
[[(0, 1/3, 2/3), (1/3, 2/3)], [(1, 0, 0), (2/3, 1/3)], [(1, 0, 0), (1, 0)]]

We see that the second equilibrium has supports of different sizes. In fact, if the first player did play \((1,0,0)\) (in other words just play the first row) the second player could play any mixture of strategies as a best response and not particularly \((2/3,1/3)\). This is because the game in consideration is now degenerate.

(Note that both of the games above are taken from Nisan et al. 2007 [pdf].)

What is a degenerate game

A bimatrix game is called nondegenerate if the number of pure best responses to a mixed strategy never exceeds the size of its support. In a degenerate game, this definition is violated, for example if there is a pure strategy that has two pure best responses (as in the example above), but it is also possible to have a mixed strategy with support size \(k\) that has \(k+1\) strategies that are a best response.

Here is an example of this:

If we consider the mixed strategy for player 2: \(y=(1/2,1/2)\), then the utility to player 1 is given by:

We see that there are 3 best responses to \(y\) and as \(y\) has support size 2 this implies that the game above is degenerate.

What does the literature say about degenerate games

The original definition of degenerate games was given in Lemke, Howson 1964 [pdf] and their definition was dependent on the labeling polytope that they used for their famous algorithm for the computation of equilibria (which is currently being implemented in Sage!). Further to this Stengel 1999 [ps] offers a nice overview of a variety of equivalent definitions.

Sadly, all of these definitions require finding a particular mixed strategy profile \((x, y)\) for which a particular condition holds. To be able to implement a test for degeneracy based on any of these definitions would require a continuous search over possible mixed strategy pairs.

In the previous example (where we take \(y=(1/2,1/2)\) we could have identified this \(y\) by looking at the utilities for each pure strategy for player 1 against \(y=(y_1, 1-y_1)\):

(\(r_i\) denotes row strategy \(i\) for player 1.) A plot of this is shown:

We can (in this instance) quickly search through values of \(y_1\) and identify the point that has the most best responses which gives the best chance of passing the degeneracy condition (\(y_1=1/2\)). This is not really practical from a generic point of view which leads to this blog post: we have identified what the particular \(x, y\) is that is sufficient to test.

A sufficient mixed strategy to test for degeneracy

The definition of degeneracy can be written as:

Def. A Normal Form Game is degenerate iff:

There exists \(x\in \Delta X\) such that \( |S(x)| < |\sigma_2| \) where \(\sigma_2\) is the support such that \( (xB)_j = \max(xB) \), for all \(j \) in \( \sigma_2\).


There exists \(y\in \Delta Y\) such that \( |S(x)| < |\sigma_1| \) where \(\sigma_1\) is the support such that \( (Ay)_i = \max(Ay) \), for all \(i \) in \( \sigma_1\).

(\(X\) and \(Y\) are the pure strategies for player 1 and 2 and \(\Delta X, \Delta Y\) the corresponding mixed strategies spaces.

The result we are implementing in Sage aims to remove the need to search particular mixed strategies \(x, y\) (a continuous search) and replace that by a search over supports (a discrete search).

Theorem. A Normal Form Game is degenerate iff:

There exists \( \sigma_1 \subseteq X \) and \( \sigma_2 \subseteq Y \) such that \( |\sigma_1| < |\sigma_2| \) and \( S(x^*) = \sigma_1 \) where \( x^* \) is a solution of \( (xB)_j = \max(xB) \), for all \(j \) in \( \sigma_2 \) (note that a valid \(x^*\) is understood to be a mixed strategy vector).


There exists \( \sigma_1 \subseteq X \) and \( \sigma_2 \subseteq Y \) such that \( |\sigma_1| > |\sigma_2| \) and \( S(y^*) = \sigma_2 \) where \( y^* \) is a solution of \( (Ay)_i = \max(Ay) \), for all \(i \) in \( \sigma_1 \).

Using the definition given above the proof is relatively straightforward but we will include it below (mainly to try and convince ourselves that we haven’t made a mistake).

We will only consider the first part of each condition (the ones for the first player). The result follows in the same way for the second player.

Proof \(\Leftarrow\)

Assume a game defined by \(A, B\) is degenerate, by the above definition without loss of generality this implies that there exists an \(x\in \Delta X\) such that \( |S(x)| < |\sigma_2| \) where \(\sigma_2\) is the support such that \( (xB)_j = \max(xB) \), for all \(j \) in \( \sigma_2\).

If we denote \(S(x)\) by \(\sigma_1\) then the definition implies that \(|\sigma_1| < |\sigma_2| \) and further more that \( (xB)_j = \max(xB) \), for all \(j \) in \( \sigma_2 \) as required.

Proof \(\Rightarrow\)

If we now assume that we have \(\sigma_1, \sigma_2, x^*\) as per the first part of the theorem then we have \(|\sigma_1|<|\sigma_2|\) and taking \(x=x^*\) implies that \(|S(x)|<|\sigma_2|\). Furthermore as \(x^*\) is a solution of \( (xB)_j = \max(xB) \) the result follows (by the definition given above).


This result implies that we simply need to consider all potential pairs of supports. Depending on the relative size of the supports we can use one of the two conditions of the result. If we ordered the supports by size the situation for the two player game looks somewhat like this:

Note that for an \(m\times n\) game there are \((2^m-1)\) potential supports for player 1 (the size of the powerset of strategy set without the empty set) and \((2^n-1)\) potential supports of for player 2. Thus the rectangle drawn above has dimension \((2^m-1)\times(2^n-1)\). Needless to say that our implementation will not be efficient (testing degeneracy is after all an NP complete problem in linear programming (see Chandrasekaran 1982 - [pdf]) but at least we have identified exactly which mixed strategy we need to test for each support pair.


  • Chandrasekaran, R., Santosh N. Kabadi, and Katta G. Murthy. “Some NP-complete problems in linear programming.” Operations Research Letters 1.3 (1982): 101-104. [pdf]
  • Lemke, Carlton E., and Joseph T. Howson, Jr. “Equilibrium points of bimatrix games.” Journal of the Society for Industrial & Applied Mathematics 12.2 (1964): 413-423. [pdf]
  • N Nisan, T Roughgarden, E Tardos, VV Vazirani Vol. 1. Cambridge: Cambridge University Press, 2007. [pdf]
  • von Stengel, B. “Computing equilibria for two person games.” Technical report. [ps]

June 25, 2015 12:00 AM

June 14, 2015

Vince Knight

Python, natural language processing and predicting funny

Every year there is a big festival in Edinburgh called the fringe festival. I blogged about this a while ago, in that post I did a very basic bit of natural language processing aiming to try and identify what made things funny. In this blog post I’m going to push that a bit further by building a classification model that aims to predict if a joke is funny or not. (tldr: I don’t really succeed but but that’s mainly because I have very little data - having more data would not necessarily guarantee success either but the code and approach is what’s worth taking from this post… 😪).

If you want to skip the brief description and go straight to look at the code you can find the ipython notebook on github here and on cloud.sagemath here.

The data comes from a series of BBC articles which reports (more or less every year since 2011?) the top ten jokes at the fringe festival. This does in fact only give 60 odd jokes to work with…

Here is the latest winner (by Tim Vine):

I decided to sell my Hoover… well it was just collecting dust.

After cleaning it up slightly I’ve thrown that all in a json file here. So in order to import the data in to a panda data frame I just run:

import pandas
df = pandas.read_json('jokes.json') # Loading the json file

Pandas is great, I’ve been used to creating my own bespoke classes for handling data but in general just using pandas does the exact right job. At this point I basically follow along with this post on sentiment analysis of twitter which makes use of the ridiculously powerful nltk library.

We can use the nltk library to ‘tokenise’ and get rid of common words:

commonwords = [e.upper() for e in set(nltk.corpus.stopwords.words('english'))] # <- Need to download the corpus: import nltk; nltk.download()
commonwords.extend(['M', 'VE'])  # Adding a couple of things that need to be removed
tokenizer = nltk.tokenize.RegexpTokenizer(r'\w+')  # To be able to strip out unwanted things in strings
string_to_list = lambda x: [el.upper() for el in tokenizer.tokenize(x) if el.upper() not in commonwords]
df['Joke'] = df['Raw_joke'].apply(string_to_list)

Note that this requires downloading one of the awesome corpuses (thats apparently the right way to say that) from nltk.

Here is how this looks:

joke = 'I decided to sell my Hoover... well it was just collecting dust.'

which gives:


We can now get started on building a classifier

Here is the general idea of what will be happening:

First of all we need to build up the ‘features’ of each joke, in other words pull the words out in to a nice easy format.

To do that we need to find all the words from our training data set, another way of describing this is that we need to build up our dictionary:

df['Year'] = df['Year'].apply(int)

def get_all_words(dataframe):
    A function that gets all the words from the Joke column in a given dataframe
    all_words = []
    for jk in dataframe['Joke']:
    return all_words

all_words = get_all_words(df[df['Year'] <= 2013])  # This uses all jokes before 2013 as our training data set.

We then build something that will tell us for each joke which of the overall words is in it:

def extract_features(joke, all_words):
    words = set(joke)
    features = {}
    for word in words:
        features['contains(%s)' % word] = (word in all_words)
    return features

Once we have done that, we just need to decide what we will call a funny joke. For this purpose We’ll use a funny_threshold and any joke that ranks above the funny_threshold in any given year will be considered funny:

funny_threshold = 5
df['Rank'] = df['Rank'].apply(int)
df['Funny'] = df['Rank'] <= funny_threshold

Now we just need to create a tuple for each joke that puts the features mentioned earlier and a classification (if the joke was funny or not) together:

df['Labeled_Feature'] = zip(df['Features'],df['Funny'])

We can now (in one line of code!!!!) create a classifier:

classifier = nltk.NaiveBayesClassifier.train(df[df['Year'] <= 2013]['Labeled_Feature'])

This classifier will take into account all the words in a given joke and spit out if it’s funny or not. It can also give us some indication as to what makes a joke funny or not:


Here is the output of that:

Most Informative Features
     contains(GOT) = True   False : True   =  2.4 : 1.0
    contains(KNOW) = True    True : False  =  1.7 : 1.0
  contains(PEOPLE) = True   False : True   =  1.7 : 1.0
     contains(SEX) = True   False : True   =  1.7 : 1.0
   contains(NEVER) = True   False : True   =  1.7 : 1.0
      contains(RE) = True    True : False  =  1.6 : 1.0
  contains(FRIEND) = True    True : False  =  1.6 : 1.0
     contains(SAY) = True    True : False  =  1.6 : 1.0
  contains(BOUGHT) = True    True : False  =  1.6 : 1.0
     contains(ONE) = True    True : False  =  1.5 : 1.0

This immediately gives us some information:

  • If your joke is about SEX is it more likely to not be funny.
  • If your joke is about FRIENDs is it more likely to be funny.

That’s all very nice but we can now (theoretically - again, I really don’t have enough data for this) start using the mathematical model to tell you if something is funny:

joke = 'Why was 10 afraid of 7? Because 7 8 9'
classifier.classify(extract_features(string_to_list(joke), get_all_words(df[df['Year'] <= 2013])))

That joke is apparently funny (the output of above is True). The following joke however is apparently not (the output of below if False):

joke = 'Your mother is ...'
print classifier.classify(extract_features(string_to_list(joke), get_all_words(df[df['Year'] <= 2013])))

As you can see in the ipython notebook it is then very easy to measure how good the predictions are (I used the data from years before 2013 to predict 2014).


Here is a plot of the accuracy of the classifier for changing values of funny_threshold:

You’ll notice a couple of things:

  • When the threshold is 0 or 1: the classifier works perfectly. This makes sense: all the jokes are either funny or not so it’s very easy for the classifier to do well.
  • There seems to be a couple of regions where the classifier does particularly poorly: just after a value of 4. Indeed there are points where the classifier does worse than flipping a coin.
  • At a value of 4, the classifier does particularly well!

Now, one final thing I’ll take a look at is what happens if I start randomly selecting a portion of the entire data set to be the training set:

Below are 10 plots that correspond to 50 repetitions of the above where I randomly sample a ratio of the data set to be the training set:

Finally (although it’s really not helpful), here are all of those on a single plot:

First of all: all those plots are basically one line of seaborn code which is ridiculously cool. Seaborn is basically magic:

sns.tsplot(data, steps)

Second of all, it looks like the lower bound of the classifiers is around .5. Most of them start of at .5, in other words they are as good as flipping a coin before we let them learn from anything, which makes sense. Finally it seems that the threshold of 4 classifier seems to be the only one that gradually improves as more data is given to it. That’s perhaps indicating that something interesting is happening there but that investigation would be for another day.

All of the conclusions about the actual data should certainly not be taken seriously: I simply do not have enough data. But, the overall process and code is what is worth taking away. It’s pretty neat that the variety of awesome python libraries lets you do this sort of thing more or less out of the box.

Please do take a look at this github repository but I’ve also just put the notebook on cloud.sagemath so assuming you pip install the libraries and get the data etc you can play around with this right in your browser:

Here is the notebook on cloud.sagemath.

June 14, 2015 12:00 AM

June 09, 2015

Michele Borassi

Performance Comparison of Different Graph Libraries

As promised in the last post, I have compared the performances of several graph libraries, in order to choose which ones should be deployed with Sagemath. Here, I provide the main results of this analysis, while more details are available on my website (see also the links below).
The libraries chosen are the most famous graph libraries written in Python, C, or C++ (I have chosen these languages because they are easier to integrate in Sagemath, using Cython). Furthermore, I have excluded NetworkX, which is already deployed with Sagemath.
First of all, I have to enforce that no graph library comparison can be completely fair, and also this comparison can be criticized, due to the large amount of available routines, to the constant evolution of libraries, and to many small differences in the outputs (for instance, one library might compute the value of a maximum s-t flow, another library might actually compute the flow, and a third one might compute all maximum flows). Despite this, I have tried to be as fair as possible, through a deeper and more detailed analysis than previous comparisons (https://graph-tool.skewed.de/performance, http://www.programmershare.com/3210372/, http://arxiv.org/pdf/1403.3005.pdf).
The first comparison deals with the number of algorithms implemented. I have chosen a set of 107 possible algorithms, trying to cover all possible tasks that a graph library should perform (avoiding easy tasks that are common to all libraries, like outputting the number of nodes, the number of edges, the neighbors of a node, etc). In some cases, two tasks were collapsed in one, if the algorithms solving these tasks are very similar (for instance, computing a maximum flow and computing a minimum cut, computing vertex betweenness and edge betweenness, etc).
The number of routines available for each library is plotted in the following chart, and a table containing all features is available in HTML or as a Google Sheet.

The results show that Sagemath has more routines than all competitors (66), closely followed by igraph (62). All other libraries are very close to each other, having about 30 routines each. Furthermore, Sagemath could be improved in the fields of neighbor similarity measures (assortativity, bibcoupling, cocitation, etc), community detection, and random graph generators. For instance, igraph contains 29 routines that are not available in Sagemath.

The second comparison analyzes the running-time of some of the algorithms implemented in the libraries. In particular, I have chosen 8 of the most common tasks in graph analysis: computing the diameter, computing the maximum flow between two vertices, finding connected components and strongly connected components, computing betweenness centrality, computing the clustering coefficient, computing the clique number, and generating a graph with the preferential attachment model. I have run each of these algorithms on 3 inputs, and I have considered the total execution time (excluding the time needed to load the graph). More details on this experiment are available here, and the results are also available in a Google Sheet.
In order to make the results more readable, I have plotted the ratio between the time needed by a given library and the minimum time needed by any library. If an algorithm was not implemented, or it needed more than 3 hours to complete, the corresponding bar is not shown.

Overall, the results show that NetworKit is the fastest library, or one of the fastest, in all routines that are implemented (apart from the generation of preferential attachment graphs, where it is very slow). Boost graph library is very close to NetworKit, and it also contains more routines. Also Sagemath is quite efficient in all tasks, apart from the computation of strongly connected components and the generation of a preferential attachment graph, where it needed more than 3 hours. However, in the latter case, the main problem was not speed but memory consumption.

In conclusion, Sagemath can highly benefit from the possibility of using algorithms from other libraries. First of all, it might improve the number of algorithms offered, especially by including igraph, and it also might improve its performance, by including Boost, NetworKit, or other fast graph libraries.

by Michele Borassi (noreply@blogger.com) at June 09, 2015 07:30 AM

June 04, 2015

Michele Borassi

Comparison of Graph Libraries

Many times, people asked me "Which is the best available graph library?", or "Which graph library should I use to compute this, or that?".
Well, personally I love to use Sage, but there are also several good alternatives. Then, the question becomes "How could we improve Sage, so that people will choose it?".

In my opinion, graph libraries are compared according to the following parameters:
  1. simplicity and documentation: people have little time, and the faster they learn how to use the library, the better;
  2. number of routines available;
  3. speed: sometimes, the input is very big, and the algorithms take much time to finish, so that a fast implementation is fundamental.
While it is very difficult to measure the first point, the others can be compared and improved. For this reason, in order to outperform other libraries, we should implement new features, and improve existing ones. You don't say!

However, this answer is not satisfactory: in principle, we could add all features available in other libraries, but this is a huge translational work, and while we are doing this work the other libraries will change, making this effort a never-ending story.

My project proposes an alternative: cooperating instead of competing. I will try to interface Sage with other libraries, and to use their algorithms when the Sage counterpart is not available, or less efficient. This way, with an affordable amount of work, we will be able to run all algorithms available in the best graph libraries!

As a first step, I have compared all the most famous C, C++, and Python graph libraries according to points 2 and 3, in order to choose which libraries should be included. The next posts will analyze the results of this comparison.

by Michele Borassi (noreply@blogger.com) at June 04, 2015 02:11 AM

Google Summer of Code: let's start!

This blog will follow my Google Summer of Code project, entitled Performance Improvements for the Graph Module of Sagemath. The complete project is available here, and related documents with partial results will be available on the same website.
In this first post, I would like to thank my mentor David Coudert and Nathann Cohen, who helped me a lot in writing this project and understanding how the graph module of Sagemath works.
With their help, and with the help of the Sage community, I hope it will be a useful and funny work! Let's start!

by Michele Borassi (noreply@blogger.com) at June 04, 2015 01:10 AM

May 29, 2015

Benjamin Hackl

Asymptotic Expressions: Motivation

\( \def\R{\mathbb{R}} \)So, as Google Summer of Code started on Monday, May 25th it is time to give a proper motivation for the project I have proposed. The working title of my project is (Multivariate) Asymptotic Expressions, and its overall goal is to bring asymptotic expressions to SageMath.

What are Asymptotic Expressions?

A motivating answer for this question comes from the theory of Taylor series. Assume that we have a sufficiently nice (in this case meaning smooth) function $f : \R \to \R$ that we want to approximate in a neighborhood of some point $x_0 \in \R$. Taylor’s theorem allows us to write $f(x) = T_n(x) + R_n(x)$ where

\[ T_n(x) = \sum_{j=0}^n \frac{f^{(j)}(x_0)}{j!}\cdot (x-x_0)^j = f(x_0) + f'(x_0)\cdot (x-x_0) + \cdots + \frac{f^{(n)}(x_0)}{n!}\cdot (x-x_0)^n,  \]

and $R_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!} \cdot (x-x_0)^{n+1}$, where $\xi$ lies in a neighborhood of $x_0$. Note that for $x\to x_0$, $R_n(x)$ “behaves like” $(x-x_0)^{n+1}$. In particular, we can certainly find a constant $C > 0$ such that $|R_n(x)| \leq C\cdot |x-x_0|^{n+1}$, or, in other words: for $x\to x_0$ the growth of the function $R_n(x)$ is bounded from above by the growth of $(x-x_0)^{n+1}$.

The idea of bounding the growth of a function by the growth of another function when the argument approaches some number (or $\infty$) is the central idea behind the big O notationFor function $f, g : \R \to \R$ we write $f(x) = O(g(x))$ for $x\to x_0$ if there is a constant $C > 0$ such that $|f(x)| \leq C\cdot |g(x)|$ for all $x$ in some neighborhood of $x_0$.

A case that is particularly important is the case of $x_0 = \infty$, that is if we want to compare and/or characterize the behavior of some function for $x\to\infty$, which is also called the functions asymptotic behavior. For example, consider the functions $\log x$, $x^3$ and $e^x$. All of them are growing unbounded for $x\to\infty$ — however, their asymptotic behavior differs. This can be seen by considering pairwise quotients of these functions: $\frac{x^3}{e^x} \to 0$ for $x\to\infty$, and therefore the asymptotic growth of $x^3$ can be bounded above by the growth of $e^x$, meaning $x^3 = O(e^x)$ for $x\to\infty$.

The analysis of a functions asymptotic behavior is important for many applications, for example when determining time and space complexity of algorithms in computer science, or for describing the growth of classes of combinatorial objects: take, for example, binary strings of length $2n$ that contain equally many zeros and ones. If $s_n$ denotes the number of such strings, then we have

\[ s_n = \binom{2n}{n} = \frac{4^n}{\sqrt{n\pi}} \left(1 + O\left(\frac{1}{n}\right)\right) \quad\text{ for } n\to\infty. \]

Expressions like these are asymptotic expressions. When we consider asymptotic expressions in only one variable, everything works out nicely as a total order is induced. But as soon as multiple variables are involved, we don’t have a total order any more. Consider, for example, $x^2 y$ and $xy^2$ when $x$ and $y$ approach $\infty$. These two elements cannot be compared to each other, which complicates computing with these expressions as they may contain multiple “irreducible” O-terms.

The following univariate and multivariate examples shall demonstrate how computing with such expressions looks like (all variables are assumed to go to $\infty$):

\[ x + O(x) = O(x),\quad x^2 \cdot (x + O(1)) = x^3 + O(x^2),\quad O(x^2) \cdot O(x^3) = O(x^5),  \]

\[ x y + O(x^2 y) = O(x^2y),\quad (y \log y + O(y)) (x^2 y + O(4^x \sqrt{x})) =  x^2 y^2 \log y + O(x^2 y^2) + O(4^x \sqrt{x} y \log y).   \]

Our plan is to provide an implementation based on which computations with these and more complicated expressions are possible.

Planned Structure

There are four core concepts of our implementation.

  • Asymptotic Growth Groups: These are multiplicative groups that contain growth elements like $x^2$, $\log x$, $2^x \cdot x \cdot \log x$. For starters, only univariate power growth groups will be implemented.
  • Asymptotic Term Monoids: These monoids contain asymptotic terms — in essence, these are summands of asymptotic terms. Apart from exact term monoids (growth elements with a coefficient), we will also implement O-term monoids as well as a term monoid for a deviation of O-terms. Asymptotic terms have (in addition to their group operation, multiplication) absorption as an additional operation: for example, O-terms are able to absorb all asymptotically “smaller” elements.
  • Mutable Poset: As we have mentioned above, due to the fact that multivariate asymptotic expressions do not have a total order with respect to their growth, we need a partially ordered set (“Poset”) that deals with this structure such that operations like absorbing terms can be performed efficiently. The mutable poset is the central data structure that asymptotic expressions are built upon.
  • Asymptotic Ring: This is our top-level structure which is also supposed to be the main interaction object for users. The asymptotic ring contains the asymptotic expressions, i.e. intelligently managed sums of asymptotic terms. All common operations shall be possible here. Furthermore, the interface should be intelligent enough such that admissible expressions from the symbolic ring can be directly converted into elements of the asymptotic ring.

Obviously, this “planned structure” is rather superficial. However, this is only to supplement the motivation for my project with some ideas on the implementation. I’ll go a lot more into the details of what I am currently implementing in the next few blog posts!


by behackl at May 29, 2015 01:34 AM

May 27, 2015

William Stein

Guiding principles for SageMath, Inc.

In February of this year (2015), I founded a Delaware C Corporation called "SageMath, Inc.".  This is a first stab at the guiding principles for the company.    It should help clarify the relationship between the company, the Sage project, and other projects like OpenDreamKit and Jupyter/IPython.

Company mission statement:

Make open source mathematical software ubiquitous.
This involves both creating the SageMathCloud website and supporting the development and distribution of the SageMath and other software, including Jupyter, Octave, Scilab, etc. Anything open source.

Company principles:

  • Absolutely all company funded software must be open source, under a GPLv3 compatible license. We are a 100% open source company.
  • Company independence and self-determination is far more important than money. A core principle is that SMI is not for sale at any price, and will not participate in any partnership (for cost) that would restrict our freedom. This means:
    • reject any offers from corp development from big companies to purchase or partner,
    • do not take any investment money unless absolutely necessary, and then only from the highest quality investors
    • do not take venture capital ever
  • Be as open as possible about everything involving the company. What should not be open (since it is dangerous):
    • security issues, passwords
    • finances (which could attract trolls)
    • private user data
What should be open:
  • aggregate usage data, e.g., number of users.
  • aggregate data that could help other open source projects improve their development, e.g., common problems we observe with Jupyter notebooks should be provided to their team.
  • guiding principles

Business model

  • SageMathCloud is freemium with the expectation that 2-5% of users pay.
  • Target audience: all potential users of cloud-based math-related software.

SageMathCloud mission

Make it as easy as possible to use open source mathematical software in the cloud.
This means:
  • Minimize onboard friction, so in less than 1 minute, you can create an account and be using Sage or Jupyter or LaTeX. Morever, the UI should be simple and streamlined specifically for the tasks, while still having deep functionality to support expert users. Also, everything persists and can be sorted, searched, used later, etc.
  • Minimize support friction, so one click from within SMC leads to a support forum, an easy way for admins to directly help, etc. This is not at all implemented yet. Also, a support marketplace where experts get paid to help non-experts (tutoring, etc.).
  • Minimize teaching friction, so everything involving software related to teaching a course is as easy as possible, including managing a list of students, distributing and collecting homework, and automated grading and feedback.
  • Minimize pay friction, sign up for a $7 monthly membership, then simple clear pay-as-you-go functionality if you need more power.

by William Stein (noreply@blogger.com) at May 27, 2015 01:03 PM

May 23, 2015

Benjamin Hackl

Google Summer of Code — Countdown

Today I received the welcome package for attending this year’s “Google Summer of Code”! Actually, it’s pretty cool; the following things were included:

  • a blue notebook with a monochromatic GSoC 15 logo (in dark blue) printed on it
  • a sticker with a colored GSoC 15 logo
  • a pen that is both a blue ballpoint pen as well as a mechanical pencil (0.5)

Here is a photo of all this stuff:



The work on our project (multivariate) Asymptotic Expressions (in cooperation with Daniel Krenn and Clemens Heuberger) begins (or rather continues) on Monday, the 25th of May. Over the course of next week (probably in a $\varepsilon$-neighborhood of Monday) I will blog about the status quo, as well as about the motivation for the project. :-)

by behackl at May 23, 2015 01:58 AM

May 04, 2015

Vince Knight

Code on cake, poker and a number theory classification web app

I have just finished writing feedback and obtaining marks for my first year students’ presentations. These presentations follow 11 weeks during which students formed companies and worked together to come up with a ‘product’ which had to involve mathematics and code (this semester comes just after 11 weeks of learning Python and Sage). In this post I’ll briefly describe some of the great things that the students came up with.

I must say that I was blown away by the standard this year. Last year the students did exceptionally well but this year the standard was even higher, I am so grateful for the effort put in by more or less everyone.

Some of the great projects included:

  • A website that used a fitted utility function (obtained from questioning family, friends, flatmates) to rank parking lots in terms of price and distance from a given venue (the website was written in Django and the function fitted using Sage).

  • A commando training app, with an actual reservist marine who is a student of ours:

  • A story based game with an original storyline stemming from the zodiac. The presentation culminated in Geraint, Jason and I (who were the audience) retaliating to their Nerf gun attack with our (hidden under the desk) Nerf guns (we had a hunch that this group would ambush us…). The game mechanics itself was coded in pure Python and the UI was almost written in Django (that was the goal but they didn’t have the time to fully implement it).

  • A Django site that had a graphical timeline of mathematics (on click you had access to a quizz and info etc…). This was one I was particularly excited about as it’s a tool I would love to use.

  • An outreach/educational package based around cryptography. They coded a variety of cyphers in Python and also put together an excellent set of teaching resources with really well drawn characters etc… They even threw in my dog Auraya (the likeness of the drawing is pretty awesome :)):

  • I ask my students to find an original way of showcasing their code. I don’t actually know the right answer to that ‘challenge’. Most students showcase the website and/or app, some will talk me through some code but this year one group did something quite frankly awesome: code on cake. Here’s some of the code they wrote for their phone app (written with kivy):

  • One group built a fully functioning and hosted web app (after taking a look at Django they decided that Flask was the way to go for this particular tool). Their app takes in a natural number and classifies it against a number of categories, go ahead and try it right now: Categorising Numbers

  • One of the more fun presentations was for a poker simulation app that uses a prime number representation of a hand of poker to simulate all possible outcomes of a given state. This work remarkably fast and immediately spits out (with neat graphics of the cards) the probability of winning given the current cards. As well as an impressive app the students presented it very well and invited me to play a game of poker (I lost, their mark was not affected…):

    Here are a couple of screen shots of the app itself:

    Home screen:

    The input card screen:

I am missing out a bunch of great projects (including an impressive actual business that I will be delighted to talk about more when appropriate). I am very grateful to the efforts put in by all the students and wish them well during their exams.

May 04, 2015 12:00 AM

April 06, 2015

Vince Knight

My 5 reasons why jekyll + github is a terrible teaching tool.

For the past year or so I have been using jekyll for all my courses. If you do not know, in a nutshell, jekyll is a ruby framework that lets you write templates for pages and build nice websites using static markdown files for your content. Here I will describe what I think of jekyll from a pedagogic point of view, in 5 main points.

1. Jekyll is terrible because the tutorial is too well written and easy to follow.

First of all, as an academic I enjoy when things are difficult to read and follow. The Jekyll tutorial can get you up and running with a jekyll site in less than 5 minutes. It is far too clear and easy to follow. This sort of clear and to the point explanation is very dangerous from a pedagogic point of view as students might stumble upon it and raise their expectations of the educational process they are going through.

In all seriousness, the tutorial is well written and clear, with a basic knowledge of the command line you can modify the base site and have a website deployed in less than 10 minutes.

2. Jekyll is terrible because it works too seamlessly with github.

First of all gh-pages takes care of the hosting. Not having to use a complicated server saves far too much time. As academics we have too much free time already, I do not like getting bored.

Github promotes the sharing and openness of code, resources and processes. Using a jekyll site in conjunction with github means that others can easily see and comment on all the materials as well as potentially improve them. This openness is dangerous as it ensures that courses are living and breathing things as opposed to a set of notes/problem sheets that sit safely in a drawer somewhere.

The fact that jekyll uses markdown is also a problem. On github anyone can easily read and send a pull request (which improves things) without really knowing markdown (let alone git). This is very terrible indeed, here for example is a pull request sent to me by a student. The student in question found a mistake in a question sheet and asked me about it, right there in the lab I just said ‘go ahead and fix it :)’ (and they did). Involving students in the process of fixing/improving their course materials has the potential for utter chaos. Furthermore normalising mistakes is another big problem: all students should be terrified of making a mistake and/or trying things.

Finally, having a personal site as a github project gives you a site at the following url:


By simply having a gh-pages branch for each class site, this will automatically be served at:


This is far too sensible and flexible. Furthermore the promotion of decentralisation of content is dangerous. If one of my class sites breaks: none of my others will be affected!!! How can I expect any free time with such a robust system? This is dangerously efficient.

3. Jekyll is terrible because it is too flexible.

You can (if you want to) include:

  • A disqus.com board to a template for a page which means that students can easily comment and talk to you about materials. Furthermore you can also use this to add things to your materials in a discussion based way, for example I have been able to far too easily to add a picture of a whiteboard explaining something students have asked.

  • Mathjax. With some escaping this works out of the box. Being able to include nicely rendered mathematics misaligns students’ expectations as to what is on the web.

  • Sage cells can be easily popped in to worksheets allowing students to immediately use code to illustrate/explain a concept.

and various others: you can just include any html/javascript etc…

This promotion of interactive and modern resources by Jekyll is truly terrible as it gets students away from what teaching materials should really be about: dusty notes in the bottom of a drawer (worked fine for me).

The flexibility of Jekyll is also really terrible as it makes me forget the restrictions imposed on me by whatever VLE we are supposed to use. This is making me weak and soft, when someone takes the choice away from me and I am forced to use the VLE, I most probably won’t be ready.

(A jekyll + github setup also implis that a wiki immediately exists for a page and I am also experimenting with a gitter.im room for each class).

4. Jekyll is terrible because it gives a responsive site out of the box.

Students should consume their materials exactly when and how we want them to. The base jekyll site cames with a basic responsive framework, here is a photo of one of my class sheets (which also again shows the disgustingly beautifully rendered mathematics):

This responsive framework works right out of the box (you can also obviously use further frameworks if you want to, see my point about flexibility) from the tutorial and this encourages students to have access to the materials on whatever platform they want whenever they want. This cannot be a good thing.

5. Jekyll is terrible because it saves me too much time.

The main point that is truly worrying about jekyll is how much time it saves me. I have mentioned this before, as academics we need to constantly make sure we do not get bored. Jekyll does not help with this.

I can edit my files using whatever system I want (I can even do this on github directly if I wanted to), I push and the website is up to date.

In the past I would have a lot of time taken up by compiling a LaTeX document and uploading to our VLE. I would sit back and worry about being bored before realising (thankfully) that I had a typo and so needed to write, delete and upload again.

Furthermore, I can easily use the github issue tracker to keep on top of to do lists etc… (which I am actually beginning to do for more or less every aspect of my life). TAs can also easily fix/improve minor things without asking me to upload whatever it is they wrote.

Github + Jekyll works seamlessly and ensures that I have more time to respond to student queries and think. This time for reflection on teaching practice is dangerous: I might choose to do things differently than how they have been done for the past 100 years.

(In case my tone is unclear: I am such a huge jekyll fan and think it is a brilliant pedagogic tool. There might well be various other static site generators and other options so please do comment about them below :))

April 06, 2015 12:00 AM

March 25, 2015

Vince Knight

A one week flipped learning environment to introduce Object Oriented Programming

This post describes a teaching activity that is run for the Cardiff MSc. programmes. The activity is revolves around a two day hackathon that gets students to use Python and object oriented programming to solve a challenge. The activity is placed within a flipped learning environment and makes use of what I feel is a very nice form of assessment (we just get to know the students).

This year is the third installment of this exercise which came as a result of the MSc advisory board requesting that object oriented programming be introduced to our MSc.

Before describing the activity itself let me just put this simple diagram that describes the flipped learning environment here (if you would like more info about it be sure to talk to Robert Talbert who has always been very helpful to me):

Description of what happens

After 3 iterations and a number of discussions about the format with Paul Harper (the director of the MSc) I think the last iteration is pretty spot on and it goes something like this:

Monday: Transfer of content

We give a brief overview of Python (you can see the slides here) up until and including basic syntax for classes.

Tuesday + Wednesday: Nothing

Students can, if they want to, read up about Python, look through videos at the website and elsewhere, look through past challenges etc… This is in effect when the knowledge transfer happens

Thursday: Flying solo followed by feedback

Students are handed a challenge of some sort (you can see the past two here). Students work in groups of 4 at attempting to solve the problem. On this day, the two postgrads (Jason and Geraint) and myself observe the groups. When we are asked questions we in general ask questions back. This sometimes leads to a fair bit of frustration but is the difficult process that makes the rest of the process worthwhile.

Here is a photo of some of the groups getting to work:

At the very end of the day (starting at 1600 for about 30 minutes with each group). During this feedback session go through the code written by each group in detail, highlighting things they are having difficulty with and agreeing on a course of action for the next day. This is the point at which the class ‘flips’ so to speak: transfer of content is done and difficulties are identified and conceptualised.

Here is a photo of Jason, Geraint and I at the end of a very long day after the feedback sessions:

The other point of this day is that we start our continuous assessment: taking notes and discussing how each group is doing:

  • Where are they progress wise?
  • What difficulties do we need to look out for?
  • How are the groups approaching the problem and working together.

Here you can see a photo of Jason in front of the board that we fill up over the 2 days with notes and comments:

Friday: Sprint finish with more assistance

On the second/last day students are given slightly more assistance from Jason, Geraint and I but are still very much left to continue with their hard work. The main difference being that when students ask questions we sometimes answer them.

Here is one group who managed to crack something quite difficult on the second day:

The final part of this day is to round all the students together and announce the marks, which brings us nicely to the assessment part of this activity.


I really enjoy assessing this activity. This is not something I say about assessment very often, but we are continuously assessing the students and are able to gain a true idea of how they do. The final piece of code is not what everything is marked on as it is in essence not terribly important.

Here is a photo of the team who did the best this year:

If I could sit with students over the 11 week period of the other courses I teach and get to know them and assess them that way, that is indeed how I would do it.


Here is a summary of how I feel this activity fits in the original diagram I had:

As you can see despite ‘being in contact’ with students for most of Thursday I would not consider this contact time in the usual sense as most of that contact is part of the assessment.

This is always a very fun (and exhausting) two days and I look forward to next year.

March 25, 2015 12:00 AM

March 24, 2015

Vince Knight

Marrying toys and students

In class yesterday we took a look at matching games. These are sometimes referred to as stable marriage problems. To have some data for us to play with I asked for some volunteers to marry. Sadly I apparently am not allowed to ask students to rank each other in class and I also do not have the authority to marry. So, like last year I used some of my office toys and asked students to rank them.

I brought three toys to class:

  • The best ninja turtle: Donatello
  • A tech deck
  • A foam football

I asked 3 students to come down and rank them and in turn I let the toys rank the students.

We discussed possible matchings with some great questions such as:

“Are we trying to make everyone as happy as possible?”

The answer to that is: no. We are simply trying to ensure that no one has an incentive to deviate from their current matching by breaking their match for someone they prefer and who also prefers them.

Here is the stable matching we found together:

Note that we can run the Gale-Shapley value using Sage:

The 3 students got to hold on to the toys for the hour and I was half expecting the football to be thrown around but sadly that did not happen. Perhaps next year.

March 24, 2015 12:00 AM

March 23, 2015

Vince Knight

Cooperative basketball in class

Today in class we repeated the game we played last year. 3 teams of 3 students took part this year and here is a photo of the aftermath:

As a class we watched the three teams attempt to free-throw as many crumpled up pieces of paper in to the bin as possible.

Based on the total number we then tried to come up with how many each subset/coalition of players would have gotten in. So for example, 2 out of 3 of the teams had one student crumple paper while the other 2 took shots. So whilst that individual did not get any in, they contributed an important part to the team effort.

Here are the characteristic functions that show what each team did:

Here is some Sage code that gives the Shapley value for each game (take a look at my class notes or at last years post to see how to calculate this):

Let us define the first game:

If you click Evaluate above you see that the Shapley value is given by:

(This one we calculated in class)

By changing the numbers above we get the following for the other two games.

  • Game 2:

  • Game 3:

This was a bit of fun and most importantly from a class point of view gave us some nice numbers to work from and calculate the Shapley value together.

If anyone would like to read about the Shapley value a bit more take a look at the Sage documentation which not only shows how to calculate it using Sage but also goes over some of the mathematics (including another formulation).

March 23, 2015 12:00 AM

March 20, 2015

Liang Ze

Character Theory Basics

This post illustrates some of SageMath’s character theory functionality, as well as some basic results about characters of finite groups.

Basic Definitions and Properties

Given a representation $(V,\rho)$ of a group $G$, its character is a map $ \chi: G \to \mathbb{C}$ that returns the trace of the matrices given by $\rho$:

A character $\chi$ is irreducible if the corresponding $(V,\rho)$ is irreducible.

Despite the simplicity of the definition, the (irreducible) characters of a group contain a surprising amount of information about the group. Some big theorems in group theory depend heavily on character theory.

Let’s calculate the character of the permutation representation of $D_4$. For each $g \in G$, we’ll display the pairs:

(The Sage cells in this post are linked, so things may not work if you don’t execute them in order.)

Many of the following properties of characters can be deduced from properties of the trace:

  1. The dimension of a character is the dimension of $V$ in $(V,\rho)$. Since $\rho(\text{Id})$ is always the identity matrix, the dimension of $\chi$ is $\chi(\text{Id})$.
  2. Because the trace is invariant under similarity transformations, $\chi(hgh^{-1}) = \chi(g)$ for all $g,h \in G$. So characters are constant on conjugacy classes, and are thus class functions.
  3. Let $\chi_V$ denote the character of $(V,\rho)$. Recalling the definitions of direct sums and tensor products, we see that

The Character Table

Let’s ignore the representation $\rho$ for now, and just look at the character $\chi$:

This is succinct, but we can make it even shorter. From point 2 above, $\chi$ is constant on conjugacy classes of $G$, so we don’t lose any information by just looking at the values of $\chi$ on each conjugacy class:

Even shorter, let’s just display the values of $\chi$:

This single row of numbers represents the character of one representation of $G$. If we knew all the irreducible representations of $G$ and their corresponding characters, we could form a table with one row for each character. This is called the character table of $G$.

Remember how we had to define our representations by hand, one by one? We don’t have to do that for characters, because SageMath has the character tables of small groups built-in:

This just goes to show how important the character of a group is. We can also access individual characters as a functions. Let’s say we want the last one:

Notice that the character we were playing with, $[4,2,0,0,0]$, is not in the table. This is because its representation $\rho$ is not irreducible. At the end of the post on decomposing representations, we saw that $\rho$ splits into two $1$-dimensional irreducible representations and one $2$-dimensional one. It’s not hard to see that the character of $\rho$ is the sum of rows 1,4 and 5 in our character table:

Just as we could decompose every representation of $G$ into a sum of irreducible representations, we can express any character as a sum of irreducible characters.

The next post discusses how to do this easily, by making use of the Schur orthogonality relations. These are really cool relations among the rows and columns of the character table. Apart from decomposing representations into irreducibles, we’ll also be able to prove that the character table is always square!

March 20, 2015 12:00 AM

March 19, 2015

Vince Knight

Playing stochastic games in class

The final blog post I am late in writing is about the Stochastic game we played in class last week. The particular type of game I am referring to is also called a Markov game where players play a series of Normal Form games, with the next game being picked from a random distribution the nature of which depends on the strategy profiles. In other words the choice of the players does not only impact on the utility gained by the players but also on the probability of what the net game will be… I blogged about this last year so feel free to read about some of the details there.

The main idea is that one stage game corresponds to this normal form game (a prisoner’s dilemma):

at the other we play:

The probability distributions, of the form \((x,1-x)\) where \(x\) is the probability with which we play the first game again are given by:

and the probability distribution for the second game was \((0,1)\). In essence the second game was an absorption game and so players would try and avoid it.

To deal with the potential for the game to last for ever we played this with a discounting factor \(\delta=1/2\). Whilst that discounting factor will be interpreted as such in theory, for the purposes of playing the game in class we used that as a probability at which the game continues.

You can see a photo of this all represented on the board:

We played this as a team game and you can see the results here:

As opposed to last year no actual duel lasted more than one round: I had a coded dice to sample at each step. The first random roll of the dice was to see if the game continued based on the \(\delta\) property (this in effect ‘deals with infinity’). The second random sample was to find out which game we payed next and if we ever went to the absorption games things finished there.

The winner was team B who in fact defected after the initial cooperation in the first game (perhaps that was enough to convince other teams they would be cooperative).

After playing this, we calculated (using some basic algebra examining each potential pure equilibria) the Nash equilibria for this game and found that there were two pure equilibria: both players Cooperating and both players defecting.

March 19, 2015 12:00 AM

March 17, 2015

Vince Knight

Incomplete information games in class

Last week my class and I looked at the basics of games with incomplete information. The main idea is that players do not necessarily know where there are in an extensive form game. We repeated a game I played last year that you can read about here.

Here is a picture of the game we played (for details take a look at the post from last year):

We played a round robing where everyone played against everyone else and you can see the results in these two notebook pages that Jason kept track off:

We see that the winner was Reg, who on both occasions of being the second player: went with the coin.

To find the Nash equilibria for this game we can translate it in to normal form game by doing the following two things:

  1. Identify the strategy sets for the players
  2. Averaging of the outcome probabilities

This gives the following strategies:


The strategies for the second player correspond to a 2-vector indexed by the information sets of the second player. In other words the first letter says what to do if the coin comes up as heads and the second letter says what to do if the coin comes up as tails:

  1. \(HH\): No matter what: play heads;
  2. \(HT\): If the coin comes up as heads: play heads. If the coin comes up as tails: play tails.
  3. \(TH\): If the coin comes up as heads: play tails. If the coin comes up as tails: play heads.
  4. \(TT\): No matter what: play tails;

Once we have done that and using the above ordering we can obtain the normal form game representation:

In class we obtained the Nash equilibria for this game by realising that the third column strategy (\(TH\): always disagree with the coin) was dominated and then carrying out some indifference analysis.

Here let us just throw it at Sage (here is a video showing and explaining some of the code):

The equilibria returned confirms what we did in class: the first player can more or less randomly (with bounds on the distribution) pick heads or tails but the second player should always agree with the coin.

March 17, 2015 12:00 AM

Discussing the game theory of walking/driving in class

Last week, in game theory class we looked at pairwise contest games. To introduce this we began by looking at the particular game that one could use to model the situation of two individuals walking or driving towards each other:

The above models people walking/driving towards each other and choosing a side of the road. If they choose the same side they will not walk/drive in to each other.

I got a coupe of volunteers to simulate this and ‘walk’ towards each other having picked a side. We very quickly arrived at one of the stage Nash equilibria: both players choosing left and/or choosing right.

I wrote a blog post about this a while ago when the BBC wrote an article about social convention. You can read that here.

We went on to compute the evolutionary stability of 3 potential stable equilibria:

  1. Everyone driving on the left;
  2. Everyone driving on the right;
  3. Everyone randomly picking a side each time.

Note that the above corresponds to the three Nash equilibria of the game itself. You can see this using some Sage code immediately (here is a video I just put together showing how one can use Sage to obtain Nash equilibria) - just click on ‘Evaluate’:

We did this calculations in two ways:

  1. From first principles using the definitions of evolutionary stability (this took a while). 2 Using a clever theoretic result that in effect does the analysis for us once and for all.

Both gave us the same result: driving on a given side of the road is evolutionarily stable whereas everyone randomly picking a side is not (a nudge in any given direction would ensure people picked a side).

This kind of corresponds to the two (poorly drawn) pictures below:

To further demonstrate the instability of the ‘choose a random side’ situation here is a plot of the actual evolutionary process (here is a video that shows what is happening):

We see that if we start by walking randomly the tiniest of mutation send everyone to picking a side.

March 17, 2015 12:00 AM

March 12, 2015

Liang Ze

Animated GIFs

I really should be posting about character theory, but I got distracted making some aesthetic changes to this blog (new icon and favicon!) and creating animations like this:


I’m not putting this in a SageCell because this could take quite a while, especially if you increase the number of frames (by changing the parameters in srange), but feel free to try it out on your own copy of Sage. It saves an animated GIF that loops forever (iterations = 0) at the location specified by savefile.

For more information, checkout the Sage reference for animated plots.

March 12, 2015 12:00 AM

March 08, 2015

Vince Knight

Playing an infinitely repeated game in class

Following the iterated Prisoner’s dilemma tournament my class I and I played last week we went on to play a version of the game where we repeated things infinitely many times. This post will briefly describe what we got up to.

As you can read in the post about this activity from last year, the way we play for an infinite amount of time (that would take a while) is to apply a discounting factor \(\delta\) to the payoffs and to interpret this factor as the probability with which the game continues.

Before I go any further (and put up pictures with the team names) I need to explain something (for the readers who are not my students).

For every class I teach I insist in spending a fair while going over a mid module feedback form that is used at Cardiff University (asking students to detail 3 things they like and don’t like about the class). One student wrote (what is probably my favourite piece of feedback ever):

“Vince is a dick… but in a good way.”

Anyway, I mentioned that to the class during my feedback-feedback session and that explains one of the team names (which I found pretty amusing):

  • Orange
  • Where’s the gun
  • We don’t know
  • Vince is a good dick

Once we had the team names set up (and I stopped trying to stop laughing) I wrote some quick Python code that we could run after each iteration:

import random
continue_prob = .5

if random.random() < continue_prob:
    print 'Game continues'
    print 'Game Over'

We started off by playing with (\delta=.5) and here are the results:

You can see the various duels here:

As you can see, very little cooperation happened this way and in fact because everyone could see what everyone else was doing Orange took advantage of the last round to create a coalition and win. We also see one particular duel that cost two teams very highly (because the ‘dice rolls’ did not really help).

After this I suggest to the class that we play again but that no one got to see what was happening to the other teams (this was actually suggested to me by students the year before). We went ahead with this and used \(delta=.25\): so the game had a less chance of carrying on.

You can see the result and duels here (this had to be squeezed on to a board that could be hidden):

Based on the theory we would expect more cooperation to be likely but as you can see this did not happen.

The tie at the end was settled with a game of Rock Paper Scissors Lizard Spock which actually gave place to a rematch of the Rock Paper Scissors Lizard Spock tournament we played earlier. Except this time Laura, lost her crown :)

March 08, 2015 12:00 AM

February 26, 2015

Sébastien Labbé

Arnoux-Rauzy-Poincaré sequences

In a recent article with Valérie Berthé [BL15], we provided a multidimensional continued fraction algorithm called Arnoux-Rauzy-Poincaré (ARP) to construct, given any vector \(v\in\mathbb{R}_+^3\), an infinite word \(w\in\{1,2,3\}^\mathbb{N}\) over a three-letter alphabet such that the frequencies of letters in \(w\) exists and are equal to \(v\) and such that the number of factors (i.e. finite block of consecutive letters) of length \(n\) appearing in \(w\) is linear and less than \(\frac{5}{2}n+1\). We also conjecture that for almost all \(v\) the contructed word describes a discrete path in the positive octant staying at a bounded distance from the euclidean line of direction \(v\).

In Sage, you can construct this word using the next version of my package slabbe-0.2 (not released yet, email me to press me to finish it). The one with frequencies of letters proportionnal to \((1, e, \pi)\) is:

sage: from slabbe.mcf import algo
sage: D = algo.arp.substitutions()
sage: it = algo.arp.coding_iterator((1,e,pi))
sage: w = words.s_adic(it, repeat(1), D)
word: 1232323123233231232332312323123232312323...

The factor complexity is close to 2n+1 and the balance is often less or equal to three:

sage: w[:10000].number_of_factors(100)
sage: w[:100000].number_of_factors(1000)
sage: w[:1000].balance()
sage: w[:2000].balance()

Note that bounded distance from the euclidean line almost surely was proven in [DHS2013] for Brun algorithm, another MCF algorithm.

Other approaches: Standard model and billiard sequences

Other approaches have been proposed to construct such discrete lines.

One of them is the standard model of Eric Andres [A03]. It is also equivalent to billiard sequences in the cube. It is well known that the factor complexity of billiard sequences is quadratic \(p(n)=n^2+n+1\) [AMST94]. Experimentally, we can verify this. We first create a billiard word of some given direction:

sage: from slabbe import BilliardCube
sage: v = vector(RR, (1, e, pi))
sage: b = BilliardCube(v)
sage: b
Cubic billiard of direction (1.00000000000000, 2.71828182845905, 3.14159265358979)
sage: w = b.to_word()
sage: w
word: 3231232323123233213232321323231233232132...

We create some prefixes of \(w\) that we represent internally as char*. The creation is slow because the implementation of billiard words in my optional package is in Python and is not that efficient:

sage: p3 = Word(w[:10^3], alphabet=[1,2,3], datatype='char')
sage: p4 = Word(w[:10^4], alphabet=[1,2,3], datatype='char') # takes 3s
sage: p5 = Word(w[:10^5], alphabet=[1,2,3], datatype='char') # takes 32s
sage: p6 = Word(w[:10^6], alphabet=[1,2,3], datatype='char') # takes 5min 20s

We see below that exactly \(n^2+n+1\) factors of length \(n<20\) appears in the prefix of length 1000000 of \(w\):

sage: A = ['n'] + range(30)
sage: c3 = ['p_(w[:10^3])(n)'] + map(p3.number_of_factors, range(30))
sage: c4 = ['p_(w[:10^4])(n)'] + map(p4.number_of_factors, range(30))
sage: c5 = ['p_(w[:10^5])(n)'] + map(p5.number_of_factors, range(30)) # takes 4s
sage: c6 = ['p_(w[:10^6])(n)'] + map(p6.number_of_factors, range(30)) # takes 49s
sage: ref = ['n^2+n+1'] + [n^2+n+1 for n in range(30)]
sage: T = table(columns=[A,c3,c4,c5,c6,ref])
sage: T
  n    p_(w[:10^3])(n)   p_(w[:10^4])(n)   p_(w[:10^5])(n)   p_(w[:10^6])(n)   n^2+n+1
  0    1                 1                 1                 1                 1
  1    3                 3                 3                 3                 3
  2    7                 7                 7                 7                 7
  3    13                13                13                13                13
  4    21                21                21                21                21
  5    31                31                31                31                31
  6    43                43                43                43                43
  7    52                55                56                57                57
  8    63                69                71                73                73
  9    74                85                88                91                91
  10   87                103               107               111               111
  11   100               123               128               133               133
  12   115               145               151               157               157
  13   130               169               176               183               183
  14   144               195               203               211               211
  15   160               223               232               241               241
  16   176               253               263               273               273
  17   192               285               296               307               307
  18   208               319               331               343               343
  19   224               355               368               381               381
  20   239               392               407               421               421
  21   254               430               448               463               463
  22   268               470               491               507               507
  23   282               510               536               553               553
  24   296               552               583               601               601
  25   310               596               632               651               651
  26   324               642               683               703               703
  27   335               687               734               757               757
  28   345               734               787               813               813
  29   355               783               842               871               871

Billiard sequences generate paths that are at a bounded distance from an euclidean line. This is equivalent to say that the balance is finite. The balance is defined as the supremum value of difference of the number of apparition of a letter in two factors of the same length. For billiard sequences, the balance is 2:

sage: p3.balance()
sage: p4.balance() # takes 2min 37s

Other approaches: Melançon and Reutenauer

Melançon and Reutenauer [MR13] also suggested a method that generalizes Christoffel words in higher dimension. The construction is based on the application of two substitutions generalizing the construction of sturmian sequences. Below we compute the factor complexity and the balance of some of their words over a three-letter alphabet.

On a three-letter alphabet, the two morphisms are:

sage: L = WordMorphism('1->1,2->13,3->2')
sage: R = WordMorphism('1->13,2->2,3->3')
sage: L
WordMorphism: 1->1, 2->13, 3->2
sage: R
WordMorphism: 1->13, 2->2, 3->3

Example 1: periodic case \(LRLRLRLRLR\dots\). In this example, the factor complexity seems to be around \(p(n)=2.76n\) and the balance is at least 28:

sage: from itertools import repeat, cycle
sage: W = words.s_adic(cycle((L,R)),repeat('1'))
sage: W
word: 1213122121313121312212212131221213131213...
sage: map(W[:10000].number_of_factors, [10,20,40,80])
[27, 54, 110, 221]
sage: [27/10., 54/20., 110/40., 221/80.]
[2.70000000000000, 2.70000000000000, 2.75000000000000, 2.76250000000000]
sage: W[:1000].balance()  # takes 1.6s
sage: W[:2000].balance()  # takes 6.4s

Example 2: \(RLR^2LR^4LR^8LR^{16}LR^{32}LR^{64}LR^{128}\dots\) taken from the conclusion of their article. In this example, the factor complexity seems to be \(p(n)=3n\) and balance at least as high (=bad) as \(122\):

sage: W = words.s_adic([R,L,R,R,L,R,R,R,R,L]+[R]*8+[L]+[R]*16+[L]+[R]*32+[L]+[R]*64+[L]+[R]*128,'1')
sage: W.length()
sage: map(W.number_of_factors, [10, 20, 100, 200, 300, 1000])
[29, 57, 295, 595, 895, 2981]
sage: [29/10., 57/20., 295/100., 595/200., 895/300., 2981/1000.]
sage: W[:1000].balance()  # takes 1.6s
sage: W[:2000].balance()  # takes 6s

Example 3: some random ones. The complexity \(p(n)/n\) occillates between 2 and 3 for factors of length \(n=1000\) in prefixes of length 100000:

sage: for _ in range(10):
....:     W = words.s_adic([choice((L,R)) for _ in range(50)],'1')
....:     print W[:100000].number_of_factors(1000)/1000.

For ten randomly generated words, the balance goes from 6 to 27 which is much more than what is obtained for billiard words or by our approach:

sage: for _ in range(10):
....:     W = words.s_adic([choice((L,R)) for _ in range(50)],'1')
....:     print W[:1000].balance(), W[:2000].balance()
12 15
8 24
14 14
5 11
17 17
14 14
6 6
19 27
9 16
12 12


[BL15]V. Berthé, S. Labbé, Factor Complexity of S-adic words generated by the Arnoux-Rauzy-Poincaré Algorithm, Advances in Applied Mathematics 63 (2015) 90-130. http://dx.doi.org/10.1016/j.aam.2014.11.001
[DHS2013]Delecroix, Vincent, Tomás Hejda, and Wolfgang Steiner. “Balancedness of Arnoux-Rauzy and Brun Words.” In Combinatorics on Words, 119–31. Springer, 2013. http://link.springer.com/chapter/10.1007/978-3-642-40579-2_14.
[A03]E. Andres, Discrete linear objects in dimension n: the standard model, Graphical Models 65 (2003) 92-111.
[AMST94]P. Arnoux, C. Mauduit, I. Shiokawa, J. I. Tamura, Complexity of sequences defined by billiards in the cube, Bull. Soc. Math. France 122 (1994) 1-12.
[MR13]G. Melançon, C. Reutenauer, On a class of Lyndon words extending Christoffel words and related to a multidimensional continued fraction algorithm. J. Integer Seq. 16, No. 9, Article 13.9.7, 30 p., electronic only (2013). https://cs.uwaterloo.ca/journals/JIS/VOL16/Reutenauer/reut3.html

by Sébastien Labbé at February 26, 2015 04:22 PM

Vince Knight

This class teaches me to not trust my classmates: An iterated prisoners dilemma in class

On Monday, in class we played an iterated prisoner’s dilemma tournament. I have done this many times (both in outreach events with Paul Harper and in this class). This is always a lot of fun but none more so than last year when Paul’s son Thomas joined us. You can read about that one here.

The format of the game is as close to that of Axelrod’s original tournament as I think it can be. I split the class in to 4 teams and we create a round robin where each team plays every other team at 8 consecutive rounds of the prisoner’s dilemma:

The utilities represent ‘years in prison’ and over the 3 matches that each team will play (against every other team) the goal is to reduce the total amount of time spent in prison.

Here are some photos from the game:

Here are the scores:

We see that ‘We will take the gun’ acquired the least total score and so they won the collection of cookies etc…

(The names followed a promise from me to let the team with the coolest name have a nerf gun… Can’t say this had the wanted effect…)

At one point during the tournament, one team actually almost declared a strategy which was cool:

We will cooperate until you defect at which point we will reevaluate

This was pretty cool as I hadn’t discussed at all what a strategy means in a repeated game (ie I had not discussed the fact that a strategy in a repeated game takes count of both play histories).

Here are all the actual duels:

You’ll also notice at the end that a coalition formed and one team agreed to defect so that they could share the prize. This happens about 50% of the time when we play this game but I never cease to be amused by it. Hopefully everyone found this fun and perhaps some even already agree with a bit of feedback I received on this course last year:

‘This class teaches me to not trust my classmates’

One of the other really cool things that happened after this class was H asking for a hand to submit a strategy to my Axelrod repository. She built a strategy called ‘Once Bitten’ that performs pretty well! You can see it here (click on ‘Blame’ and you can see the code that she wrote).

(Big thanks to Jason for keeping track of the scores and to Geraint for helping and grabbing some nice pictures)

February 26, 2015 12:00 AM

February 15, 2015

Liang Ze

The Group Ring and the Regular Representation

In the previous post, we saw how to decompose a given group representation into irreducibles. But we still don’t know much about the irreducible representations of a (finite) group. What do they look like? How many are there? Infinitely many?

In this post, we’ll construct the group ring of a group. Treating this as a vector space, we get the regular representation, which turns out to contain all the irreducible representations of $G$!

The group ring $FG$

Given a (finite) group $G$ and a field $F$, we can treat each element of $G$ as a basis element of a vector space over $F$. The resulting vector space generated by $g \in G$ is

Let’s do this is Sage with the group $G = D_4$ and the field $F = \mathbb{Q}$:

(The Sage cells in this post are linked, so things may not work if you don’t execute them in order.)

We can view $v \in FG$ as vector in $F^n$, where $n$ is the size of $G$ :

Here, we’re treating each $g \in G$ as a basis element of $FG$

Vectors in $FG$ are added component-wise:

Multiplication as a linear transformation

In fact $FG$ is also a ring (called the group ring), because we can multiply vectors using the multiplication rule of the group $G$:

That wasn’t very illuminating. However, treating multiplication by $v \in FG$ as a function

one can check that each $T_v$ is a linear transformation! We can thus represent $T_v$ as a matrix whose columns are $T_v(g), g \in G$:

The regular representation

We’re especially interested in $T_g, g \in G$. These are invertible, with inverse $T_{g^{-1}}$, and their matrices are all permutation matrices, because multiplying by $g \in G$ simply permutes elements of $G$:

Define a function $\rho_{FG}$ which assigns to each $g\in G$ the corresponding $T_g$:

Then $(FG,\rho_{FG})$ is the regular representation of $G$ over $F$.

The regular representation of any non-trivial group is not irreducible. In fact, it is a direct sum of all the irreducible representations of $G$! What’s more, if $(V,\rho)$ is an irreducible representation of $G$ and $\dim V = k$, then $V$ occurs $k$ times in the direct-sum decomposition of $FG$!

Let’s apply the decomposition algorithm in the previous post to $(FG,\rho_{FG})$ (this might take a while to run):

So the regular representation of $D_4$ decomposes into four (distinct) $1$-dim representations and two (isomorphic) $2$-dim ones.

Building character

We’ve spent a lot of time working directly with representations of a group. While more concrete, the actual matrix representations themselves tend to be a little clumsy, especially when the groups in question get large.

In the next few posts, I’ll switch gears to character theory, which is a simpler but more powerful way of working with group representations.

February 15, 2015 12:00 AM

February 02, 2015

Liang Ze

Decomposing Representations

In this post, we’ll implement an algorithm for decomposing representations that Dixon published in 1970.

As a motivating example, I’ll use the permutation matrix representation of $D_4$ that we saw in an earlier post. To make the code more generally applicable, let’s call the group $G$ and the representation $\rho$:

(The Sage cells in this post are linked, so things may not work if you don’t execute them in order.)

We’ll see that this is decomposable, and find out what its irreducible components are.

Unitary representations

A short remark before we begin: The algorithm assumes that $\rho$ is a unitary representation

i.e. for all $g \in G$,

where $A*$ is the conjugate transpose of a matrix $A$. For $G$ a finite group, all representations can be made unitary under an appropriate change of basis, so we need not be too concerned about this. In any case, permutation representations are always unitary, so we can proceed with our example.

Finding non-scalar, commuting matrices

At the end of the previous post we saw that in order to decompose a representation $(V,\rho)$, it is enough to find a non-scalar matrix $T$ that commutes with $\rho(g)$ for every $g \in G$. This first step finds a Hermitian non-scalar $H$ that commutes with $\rho(G)$ (if there is one to be found).

Let $E_{rs}$ denote the $n \times n$ matrix with a $1$ in the $(r,s)$th entry and zeros everywhere else. Here $n$ is the dimension of $V$ in the representation $(V,\rho)$. Define

then the set of matrices $H_{rs}$ forms a Hermitian basis for the $n \times n$ matrices over $\mathbb{C}$.

Now for each $r,s$, compute the sum

Observe that $H$ has the following properties:

  • it is hermitian
  • it commutes with $\rho(g)$ for all $g \in G$

If $\rho$ is irreducible, then $H$ is a scalar matrix for all $r,s$. Otherwise, it turns out that there will be some $r,s$ such that $H$ is non-scalar (this is due to the fact that the $H_{rs}$ matrices form a basis of the $n \times n$ matrices$).

Let’s test this algorithm on our permutation representation of $D_4$:

We get a non-scalar $H$! So the permutation representation of $D_4$ is reducible!

Using $H$ to decompose $\rho$

Our next step is to use the eigenspaces of $H$ to decompose $\rho$. At the end of the previous post, we saw that $\rho(g)$ preserves the eigenspaces of $H$, so we need only find the eigenspaces of $H$ to decompose $\rho$.

Since $H$ is hermitian, it is diagonalizable, so its eigenvectors form a basis of $V$. We can find this basis by computing the Jordan decomposition of $H$:

Finally, we observe that $P^{-1} \rho(g) P$ has the same block-diagonal form for each $g \in G$:

We have thus decomposed $\rho$ into two 1-dimensional representations and one 2-dimensional one!

Decomposing into irreducibles

Finally, to get a decomposition into irreducibles, we can apply the algorithm recursively on each of the subrepresentations to see if they further decompose.

Here’s a stand-alone script that decomposes a representation into its irreducible components:

Getting all irreducible representations

Now we know how to test for irreducibility and decompose reducible representations. But we still don’t know how many irreducible representations a group has.

It turns out that finite groups have finitely many irreducible representations! In the next post, we’ll construct a representation for any finite group $G$ that contains all the irreducible representations of $G$.

February 02, 2015 12:00 AM

January 26, 2015

Liang Ze

Irreducible and Indecomposable Representations

Following up from the questions I asked at the end of the previous post, I’ll define (ir)reducible and (in)decomposable representations, and discuss how we might detect them. Unlike previous posts, this post will have just text, and no code. This discussion will form the basis of the algorithm in the next post.


In the previous post, I showed how to form the direct sum $(V_1 \oplus V2,\rho)$ of two representations $(V_1,\rho_1)$ and $(V_2,\rho_2)$. The matrices given by $\rho$ looked like this:

A representation $(V,\rho)$ is decomposable if there is a basis of $V$ such that each $\rho(g)$ takes this block diagonal form. If $(V,\rho)$ does not admit such a decomposition, it is indecomposable.

Equivalently, $(V,\rho)$ is decomposable if there is an invertible matrix $P$ such that for all $g\in G$,

and indecomposable otherwise. Here, $P$ is a change of basis matrix and conjugating by $P$ changes from the standard basis to the basis given by the columns of $P$.


Notice that if $\rho(g)$ were block diagonal, then writing $v \in V$ as ${v_1 \choose v_2}$, where $v_1$ and $v_2$ are vectors whose dimensions agree with the blocks of $\rho(g)$, we see that

Let $V_1$ be the subspace of $V$ corresponding to vectors of the form ${v_1 \choose 0}$, and $V_2$ be the subspace of vectors of the form ${0 \choose v_2}$. Then for all $g \in G, v \in V_i$,

Now suppose instead that for all $g \in G, \rho(g)$ has the block upper-triangular form

where $ * $ represents an arbitrary matrix (possibly different for each $g \in G$). If $*$ is not the zero matrix for some $g$, then we will still have $\rho(g) v \in V_1 \,\, \forall v \in V_1$, but we no longer have $\rho(g) v \in V_2 \,\, \forall v \in V_2$. In this case, we say that $V_1$ is a subrepresentation of $V$ whereas $V_2$ is not.

Formally, if we have a subspace $W \subset V$ such that for all $g \in G, w \in W$,

then $W$ is a $G$-invariant subspace of $V$, and $(W,\rho)$ is a subrepresentation of $(V,\rho)$.

Any representation $(V,\rho)$ has at least two subrepresentations: $(0,\rho)$ and $(V,\rho)$. If there are no other subrepresentations, then $(V,\rho)$ is irreducible. Otherwise, it is reducible.

Equivalently, $(V,\rho)$ is reducible if there is an invertible matrix $P$ such that for all $g \in G$,

and irreducible otherwise.

Maschke’s Theorem

Note that a decomposable representation is also reducible, but the converse is not generally true. (Equivalently: an irreducible representation is also indecomposable, but the converse is not generally true.) Maschke’s Theorem tells us that the converse is true over fields of characteristic zero! In other words:

Suppose $V$ is a vector space over a field of characteristic zero, say $\mathbb{C}$, and $(V,\rho)$ has a subrepresentation $(W_1,\rho)$. Then there is a subspace $W_2$ (called the direct complement of $W_1$) such that $V = W_1 \oplus W_2$.

Since we will be working over $\mathbb{C}$, we can thus treat (in)decomposability as equivalent to (ir)reducibility. To understand representations of $G$, we need only understand its irreducible representations, because any other representation can be decomposed into a direct sum of irreducibles.

Schur’s Lemma

How may we detect (ir)reducible representations? We’ll make use of the following linear algebraic properties:

Given an eigenvalue $\lambda$ of a matrix $A \in \mathbb{C}^{n \times n}$, its $\lambda$-eigenspace is

Clearly, each eigenspace is an invariant subspace of $A$. If we have another matrix $B \in \mathbb{C}^{n \times n}$ such that $AB = BA$, then $B$ preserves the eigenspaces of $A$ as well. To see this, take $v \in E_\lambda$, then

so $E_\lambda$ is also an invariant subspace of $B$!

Now suppose we have a representation $(V,\rho)$ and a linear map $T:V \to V$ such that for all $g \in G, v \in V$,

Treating $T$ as a matrix, this is equivalent to saying that $\rho(g)T = T\rho(g)$ for all $g \in G$. In that case, the eigenspaces of $T$ are $G$-invariant subspaces, and will yield decompositions of $(V,\rho)$ if they are not the whole of $V$. But if $E_\lambda = V$, then $Tv = \lambda v$ for all $v \in V$, so in fact $T = \lambda I$, where $I$ is the identity matrix. We have thus shown a variant of Schur’s lemma:

If $(V,\rho)$ is irreducible, and $\rho(g) T = T \rho(g)$ for all $g \in G$, then $T =\lambda I$ for some $\lambda$.

We already know that scalar matrices (i.e. matrices of the form $\lambda I$) commute with all matrices. If $(V,\rho)$ is irreducible, this result says that there are no other matrices that commute with all $\rho(g)$. The converse is also true:

If $(V,\rho)$ is a reducible, then there is some $T \neq \lambda I$ such that $\rho(g) T = T\rho(g)$ for all $g \in G$.

I won’t prove this, but note that if $V$ has a decomposition $W_1 \oplus W_2$, then the projection onto either $W_i$ will have the desired properties. If we have such a $T$, then its eigenspaces will give a decomposition of $(V,\rho)$. This will be the subject of the next post.

January 26, 2015 12:00 AM

January 24, 2015

Liang Ze

Direct Sums and Tensor Products

In this short post, we will show two ways of combining existing representations to obtain new representations.


In the previous post, we saw two representations of $D_4$: the permutation representation, and the representation given in this Wikipedia example. Let’s first define these in Sage:

(The Sage cells in this post are linked, so things may not work if you don’t execute them in order.)

Direct Sums

If $(V_1,\rho_1), (V_2,\rho_2)$ are representations of $G$, the direct sum of these representations is $(V_1 \oplus V_2, \rho)$, where $\rho$ sends $g \in G$ to the block diagonal matrix

Here $\rho_1(g), \rho_2(g)$ and the “zeros” are all matrices.

It’s best to illustrate with an example. We can define a function direct_sum in Sage that takes two representations and returns their direct sum.

Tensor products

We can also form the tensor product $(V_1 \otimes V_2,\rho)$, where $\rho$ sends $g \in G$ to the Kronecker product of the matrices $\rho_1(g)$ and $\rho_2(g)$.

We define a function tensor_prod that takes two representations and returns their tensor product.

Observe that

  • $\dim V_1 \oplus V_2 = \dim V_1 + \dim V_2$,
  • $\dim V_1 \otimes V_2 = \dim V_1 \times \dim V_2$,

which motivates the terms direct sum and tensor product.

We can keep taking direct sums and tensor products of existing representations to obtain new ones:

Decomposing representations

Now we know how to build new representations out of old ones. One might be interested in the inverse questions:

  1. Is a given representation a direct sum of smaller representations?
  2. Is a given representation a tensor product of smaller representations?

It turns out that Q1 is a much more interesting question to ask than Q2.

A (very poor) analogy of this situation is the problem of “building up” natural numbers. We have two ways of building up new integers from old: we can either add numbers, or multiply them. Given a number $n$, it’s easy (and not very interesting) to find smaller numbers that add up to $n$. However, finding numbers whose product is $n$ is much much harder (especially for large $n$) and much more rewarding. Prime numbers also play a special role in the latter case: every positive integer has a unique factorization into primes.

The analogy is a poor one (not least because the roles of “sum” and “product” are switched!). However, it motivates the question

  • What are the analogues of “primes” for representations?

We’ll try to answer this last question and Q1 in the next few posts, and see what it means for us when working with representations in Sage.

January 24, 2015 12:00 AM

January 20, 2015

Liang Ze

Representation Theory in Sage - Basics

This is the first of a series of posts about working with group representations in Sage.

Basic Definitions

Given a group $G$, a linear representation of $G$ is a group homomorphism $\rho: G \to \mathrm{GL}(V)$.

For our purposes, we will assume that $G$ is a finite group and $V$ is an $n$-dimensional vector space over $\mathbb{C}$. Then $\mathrm{GL}(V)$ is isomorphic to the invertible $n \times n$ matrices over $\mathbb{C}$, which we will denote $\mathrm{GL}_n \mathbb{C}$.

So a representation is just a function that takes group elements and returns invertible matrices, such that

where we have matrix multiplication on the right-hand side.

Various authors refer to the map $\rho$, the vector space $V$, or the tuple $(V,\rho)$ as a representation; this shouldn’t cause any confusion, as it’s usually clear from context whether we are referring to a map or a vector space. When I need to be extra precise, I’ll use $(V,\rho)$.

Some simple examples

Trivial representation

The simplest representation is just the trivial representation that sends every element of $G$ to the identity matrix (of some fixed dimension $n$). Let’s do this for the symmetric group $S_3$:

(The Sage cells in this post are linked, so things may not work if you don’t execute them in order.)

We can verify that this is indeed a group homomorphism (warning: There are 6 elements in $S_3$, which means we have to check $6^2 = 36$ pairs!):

Permutation representation

This isn’t very interesting. However, we also know that $S_3$ is the group of permutations of the 3-element set {$1,2,3$}. We can associate to each permutation a permutation matrix. Sage already has this implemented for us, via the method matrix() for a group element g:

Qn: From the permutation matrix, can you tell which permutation $g$ corresponds to?

We can again verify that this is indeed a representation. Let’s not print out all the output; instead, we’ll only print something if it is not a representation. If nothing pops up, then we’re fine:

Defining a representation from generators

We could define permutation representations so easily only because Sage has them built in. But what if we had some other representation that we’d like to work with in Sage? Take the dihedral group $D_4$. Wikipedia tells us that this group has a certain matrix representation. How can we recreate this in Sage?

We could hard-code the relevant matrices in our function definition. However, typing all these matrices can be time-consuming, especially if the group is large.

But remember that representations are group homomorphisms. If we’ve defined $\rho(g)$ and $\rho(h)$, then we can get $\rho(gh)$ simply by multiplying the matrices $\rho(g)$ and $\rho(h)$! If we have a set of generators of a group, then we only need to define $\rho$ on these generators. Let’s do that for the generators of $D_4$:

We see that $D_4$ has a generating set of 2 elements (note: the method gens() need not return a minimal generating set, but in this case, we do get a minimal generating set). Let’s call these $r$ and $s$. We know that elements of $D_4$ can be written $r^is^j$, where $i = 0,1,2,3$ and $j = 0,1$. We first run through all such pairs $(i,j)$ to create a dictionary that tells us which group elements are given by which $(i,j)$:

Now for $g = r^i s^j \in D_4$, we can define $\rho(g) = \rho(r)^i \rho(s)^j$ and we will get a representation of $D_4$. We need only choose the matrices we want for $\rho(r)$ and $\rho(s)$.

$r$ and $s$ correspond to $R_1$ and $S_0$, resp., in the Wikipedia example, so let’s use their matrix representations to generate our representation:

One can verify that this does indeed give the same matrices as the Wikipedia example, albeit in a different order.

We can do better!

All the representations we’ve defined so far aren’t very satisfying! For the last example, we required the special property that all elements in $D_4$ have the form $r^i s^j$. In general, it isn’t always easy to express a given group element in terms of the group’s generators (this is known as the word problem).

We’ve also been constructing representations in a rather ad-hoc manner. Is there a more general way to construct representations? And how many are representations are there?

In the next post, I’ll run through two simple ways of combining existing representations to get new ones: the direct sum and the tensor product. I’ll also define irreducible representations, and state some results that will shed some light on the above questions.

January 20, 2015 12:00 AM

January 17, 2015

Liang Ze

Subgroup Explorer

To sum up the subgroup lattice series, I’ve written a subgroup lattice generator. It’s powered by Sage and GAP, and allows you to view the lattice of subgroups or subgroup conjugacy classes of a group from your browser.

Normal subgroups are colored green. Additionally, the center is blue while the commutator subgroup is pink.

Showing the full subgroup lattice can get messy for large groups. If the option Conjugacy classes is selected, the viewer only shows the conjugacy classes of subgroups (i.e. all subgroups that are conjugate are combined into a single vertex).

The edge labels indicate how many subgroups of one conjugacy class a given representative subgroup of another conjugacy class contains, or how many subgroups of one conjugacy class a given representative subgroup of another conjugacy class is contained by. The labels are omitted if these numbers are 1. The edge colors indicate whether the subgroups in the “smaller” conjugacy class are normal subgroups of those in “larger” conjugacy class.

For instance, the group C15 : C4 (of order 60; the colon stands for semi-direct product and is usually written $\rtimes$) contains 5 subgroups isomorphic to C3 : C4, which in turn contains 3 subgroups isomorphic to C4 and 1 subgroup isomorphic to C6 (the 5 belows to another edge). The edge colors indicate that C6 is a normal subgroup of C3 : C3 whereas C4 is not. For further information on group descriptors, click here.

As mentioned in the previous post, labelling the edges requires taking apart the poset plotting code. I had to extract the Hasse diagram of a poset as a graph and modifying the edge labels directly. This explains why the code is much longer than in previous posts.

Finally, while verifying the results of this program, I found an error in this book! The correction has been pencilled in. The original number printed was 1. A5 Lattice

This is the last post on subgroup lattices. My next post will be the start of a new series on doing represntation theory in Sage.

January 17, 2015 12:00 AM

December 27, 2014

Liang Ze

Lattice of Subgroups III - Coloring Edges

This post will cover the coloring of edges in the lattice of subgroups of a group.

Lattice of subgroups of $C3:C8$

Generating small groups

As we’ve done in previous posts, let’s choose a group and generate its lattice of subgroups. This can be done by referring to this list of constructions for every group of order less than 32 . These instructions allow us to construct every group on Wikipedia’s list of small groups!

For this post, we’ll use $G = C_3 \rtimes C_8$ (or $\mathbb{Z}_3 \rtimes \mathbb{Z}_8$). First, we’ll generate $G$ and display it’s poset of subgroups. For simplicity, we won’t color the vertices.

(The Sage cells in this post are linked, so things may not work if you don’t execute them in order.)

Coloring edges

In the previous post, we colored vertices according to whether the corresponding subgroup was normal (or abelian, or a Sylow subgroup, etc.) These are properties that depend only on each individual subgroup.

However, suppose we want to see the subnormal series of the group. A subnormal series is a series of subgroups where each subgroup is a normal subgroup of the next group in the series. Checking whether a particular series of subgroups is a subnormal series requires checking pairs of subgroups to see whether one is a normal subgroup of the other. This suggests that we color edges according to whether one of its endpoints is a normal subgroup of the other endpoint.

The edges of the Hasse diagram of a poset are the pairs $(h,k)$ where $h$ is covered by $k$ in the poset. This means that $h < k$, with nothing else in between. We thus obtain all the edges of a Hasse diagram by calling P.cover_relations() on the poset $P$.

To color edges of a graph, we create a dictionary edge_colors:

Up next…

This is the last post describing relatively simple things one can do to visualize subgroup lattices (or more generally, posets) in Sage. In the next post, I’ll write code to label edges, but will probably skip explaining it, as doing this requires extracting the Hasse diagram of a poset as a graph and modifying the edge labels. Also, subgroup lattices tend to get unwieldy for large groups. In the next post, we’ll restrict our attention to conjugacy classes of subgroups, rather than all subgroups.

After that, I hope to write a bit about doing some simple representation theory in Sage.

December 27, 2014 12:00 AM

December 25, 2014

Liang Ze

Holiday Harmonograph

(Guest post from the Annals of Harmonography) harmonograph

When it’s snowing outside (or maybe not),

And your feet are cold (or maybe hot),

When it’s dark as day (or bright as night),

And your heart is heavy (and head is light),

What should you do (what should you say)

To make it all right (to make it okay)?


Just pick up a pen (a pencil will do),

Set up a swing (or three, or two),

And while the world spins (or comes to a still),

In your own little room (or on top of a hill),

Let your pendulum sway (in its time, in its way),

And watch as the pen draws your worries away!



(Click inside the colored box to choose a color. Then click outside and watch it update.)

  • 7 celebrities and their harmonographs
  • What your harmonograph says about you
  • 10 tips for a happier harmonograph
  • Harmonograph secrets… revealed!

December 25, 2014 12:00 AM

November 14, 2014

William Stein

SageMathCloud Notifications are Now Better

I just made live a new notifications systems for  SageMathCloud, which I spent all week writing.  

These notifications are what you see when you click the bell in the upper right.   This new system replaces the one I made live two weeks ago.     Whenever somebody actively *edits* (using the web interface) any file in any project you collaborate on, a notification will get created or updated.    If a person *comments* on any file in any project you collaborate on (using the chat interface to the right), then not only does the notification get updated, there is also a little red counter on top of the bell and also in the title of the  SageMathCloud tab.   In particular, people will now be much more likely to see the chats you make on files.

NOTE: I have not yet enabled any sort of daily email notification summary, but that is planned. 

Some technical details:  Why did this take all week?  It's because the technology that makes it work behind the scenes is something that was fairly difficult for me to figure out how to implement.  I implemented a way to create an object that can be used simultaneously by many clients and supports realtime synchronization.... but is stored by the distributed Cassandra database instead of a file in a project.   Any changes to that object get synchronized around very quickly.   It's similar to how synchronized text editing (with several people at once) works, but I rethought differential synchronization carefully, and also figured out how to synchronize using an eventually consistent database.    This will be useful for implementing a lot other things in SageMathCloud that operate at a different level than "one single project".  For example, I plan to add functions so you can access these same "synchronized databases" from Python processes -- then you'll be able to have sage worksheets (say) running on several different projects, but all saving their data to some common synchronized place (backed by the database).   Another application will be a listing of the last 100 (say) files you've opened, with easy ways to store extra info about them.    It will also be easy to make account and project settings more realtime, so when you change something, it automatically takes effect and is also synchronized across other browser tabs you may have open.   If you're into modern Single Page App web development, this might remind you of Angular or React or Hoodie or Firebase -- what I did this week is probably kind of like some of the sync functionality of those frameworks, but I use Cassandra (instead of MongoDB, say) and differential synchronization.  

I BSD-licensed the differential synchronization code  that I wrote as part of the above. 

by William Stein (noreply@blogger.com) at November 14, 2014 02:31 PM