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 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 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 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 18, 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 18, 2015 11:01 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.
[DHS2013]Delecroix, Vincent, Tomás Hejda, and Wolfgang Steiner. “Balancedness of Arnoux-Rauzy and Brun Words.” In Combinatorics on Words, 119–31. Springer, 2013.
[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).

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 20, 2015

Vince Knight

An iterated prisoners dilemma on github

In the 1980s Axelrod ran a computer tournament inviting people to contribute code that specified strategies in an iterated prisoner’s dilemma tournament. I have just finished putting the final pieces on a Python repository on github ( to carry out the same tournament and would be delighted for people to contribute strategies via pull request (or indeed via any way possible: just get in touch). In this post I’ll describe the process of adding a strategy to the repository (the first 3 minutes of a video at the end of this post show you exactly what you are required to do).

For a good overview of the iterated prisoner’s dilemma take a look at this page about Axelrod’s tournament but in a nutshell the idea is that two players (prisoners) repeatedly play the following game:

If in a particular round they both cooperate (first row/column) they both accrue 2 years in prison. If one defects (second row/column) and the other cooperates: the defector gets 0 extra years in prison and the cooperator 5. If they both defect they each accrue 4 years in prison.

Axelrod’s tournament invited contribution of strategies that took account of the history of both players over several rounds. Thus a strategy that punished defectors would perhaps wait until a defection to do that (cooperating until then). The tournament was a round robin with the lowest total/mean years in prison being deemed the winner.

This tournament has often been used to describe how cooperation can emerge in a population: the tit for tat strategy (which starts by cooperating and then simply repeats the previous action) won! (In fact it won both times as the tournament was repeated.)

I have put together a github repository that allows anyone to contribute a strategy using Python to this tournament. You can find it here:

At present I have only implemented 6 strategies and you can see the result here:

To contribute you really only need to write very simple python code. Here is the code for the tit for tat strategy:

from axelrod import Player

class TitForTat(Player):
    A player starts by cooperating and then mimics previous move by opponent.
    def strategy(self, opponent):
        Begins by playing 'C':
        This is affected by the history of the opponent: the strategy simply repeats the last action of the opponent
            return opponent.history[-1]
        except IndexError:
            return 'C'

    def __repr__(self):
        The string method for the strategy.
        return 'Tit For Tat'

That is more or less just 3 lines of code. I’ll now briefly describe adding another strategy to this repository (and I will do it entirely using the github web interface).

I am going to add a strategy called: ‘alternator’ which simply alternates strategies.

First I navigate to this url:

Here I can just click on the + at the top of the page (after: Axelrod/axelrod/strategies/+). As this is my own github repository I can just immediately start creating a file, anyone else would be taken the github process of forking the repository:

After writing the code I simply scroll down to the bottom where I am able to commit the change but others would be able to submit a pull request:

Let us take a look at the actual code I wrote for the alternator class:

from axelrod import Player

class Alternator(Player):
    A player who alternates between cooperating and defecting
    def strategy(self, opponent):
        Alternate 'C' and 'D'
        if self.history[-1] == 'C':
            return 'D'
        return 'C'

    def __repr__(self):
        The string method for the strategy:
        return 'Alternator'

This just inherits from a Player class I have created previously and all it really requires are two methods:

  • strategy: this takes in the player itself (self) and the opponent and must return either C or D. Take a look through the other strategies to see how this can be written.
  • __repr__: this just returns what we want the strategy to look like when printed out.

After writing the code itself we also need to modify the file in the strategies directory. Here I have added the relevant lines:

from cooperator import *
from defector import *
from grudger import *
from rand import *
from titfortat import *
from gobymajority import *
from alternator import *  # <- Adding this line

strategies = [
        Alternator,  # And adding this line

Now if you want to be super awesome you can also add a test for your strategy (this helps keep things organised and working but do not let this be an obstacle to contributing a strategy). If you want to see what the test for the alternator looks like you can find it here (the README has info as to how to run the tests).

The latest results (which as of the time of writing now only include the alternator as an extra strategy - EDIT Just before publishing this Geraint sent me a pull request for another strategy: this image is the live one from the github repo so will have more and more strategies) can be found seen here:

If you clone this repository you can obtain that plot by running:

$ python

TLDR: Please contribute a strategy to this Python based version of Axelrod’s tournament on github:

This would hopefully be of interest if:

  • You are slightly interested in Axelrod’s work
  • You like Python
  • You like Game Theory
  • You have never had a pull request accepted before

Here is a short video showing exactly how to contribute using nothing else than the github web interface (the first 3 minutes are all you really need to do):

February 20, 2015 12:00 AM

February 15, 2015

Vince Knight

On the worst play in superbowl history

In my last game theory class I gave students the choice between starting to rigorously look at extensive form games or to look at the last play of this years superbowl. No one objected to looking at the superbowl play (which has been called ‘the worst play in superbowl history’) as an opportunity to go over all the normal form game theory we have seen so far.

First of all I had to give a brief overview of American football so I drew a couple of pictures on the board, here is more or less what I drew:

I highlighted that the offensive team (the blue guys in my picture) can basically at any point do one of two things:

  • Pass: the quarter back throws the ball down field. This works better when there is more room down the field for guys to get open but in general is a high risk high reward player.
  • Run: the quarter back hands the ball to the running back. This in general ensures a small gain of yards and is less risky (but has less reward).

The defence can (at a very simplistic level) set up in two ways:

  • Defend the pass: they might drop some of guys (red circles) ‘out of the box’ in a way as to cover blue guys trying to get open.
  • Defend the run: they might drop more guys (red crosses) ‘in to the box’, leaving less people to cover the pass but having more people ready to stop the running back from gaining yards.

This immediately lets us set up (in a simplistic way) every phase of play in American Football as a two player normal form game where the Offense has strategy set: \(S_O=\{\text{P},\text{R}\}\) and the Defense has strategy set: \(S_D=\{\text{DP},\text{DR}\}\).

Now, to return to what is being called ‘the worst play in superbowl history’:

  • The Seahawks have one of the best running backs in the game (funny interview of his here);
  • They had the ball very close to the score line with downs (attempts) to spare.

Everyone seems to be saying that ‘not rushing’ was an idiotic thing to do. This article by the Economist explains things pretty well (including the back story), but in my class on Friday I thought I would go through the mathematics of it all.

First of all let us build some basic utilities (to be clear I am more or less pulling these out of a hat although I will carry out some Monte Carlo simulation at the end of this blog post).

  • Let us assume (because Lynch is awesome) that on average running against a run defence would work 60% of the time;
  • Running against a pass defence would work 98% of the time;
  • Passing against a rush defence would work 90% of the time;
  • Passing against a pass defence would work 50% of the time.

Our normal form game then looks like this (assuming that the game is zero sum):

If the Seahawks were 100% going to run the ball, in such a way as that the defence were in no doubt then the obvious best response is to use a rush defence (in reality this is actually what happened and the defensive player that was fairly isolated just made an amazing play). Then, the chance of the Seahawks scoring is 60%.

What does game theory say should happen?

First of all what should the offence do to ensure that the defense sometimes defends the pass? The theory tells us that the offence should make the defense indifferent. So let us assume that the offence runs \(x\) of the time so that the player is playing mixed strategy: \(\sigma_{\text{O}}=(x,1-x)\). The Nash equilibrium strategy for the offence would then obey the following equation:

the above corresponds to:

which has solution: \(x=20/30\approx.51\). So to make the defense indifferent, the offense should run only 51% of the time!

This in itself is not a Nash equilibria: we need to calculate the strategy that the defense should play so as to make the offense indifferent. We let \(\sigma_{\text{D}}=(y,1,y)\) and now need to solve:

the above corresponds to:

which has solution: \(y=8/13\approx.62\). So to make the offense indifferent, the defense should defend the run 62% of the time.

How good is the Nash equilibria?

We discussed earlier that if the offense passed with 0 probability (so in fact always ran in that situation) then the probability of scoring would be 60% (because the defense would know and just defend the run). The Nash equilibria we have calculated ensured that the play calling is balanced: which ‘keeps the defense honest’. What is the effect of this on offenses chances of scoring?

To find this out we can simply calculate the utility at the Nash equilibria. We can do this by calculating:

or because of the indifference ensured earlier we can just calculate:

So the Nash Equilibria makes things \(74.62/60\approx 1.24\) times better for the offense. Thus, in way, by picking a strategy in that particular instance that ended up not paying off, the offense ensured a larger long term chance of scoring in a similar situation.

Some sensitivity analysis

I have completely picked numbers more or less out of a hat. Thankfully we can combine the above analysis with some Monte Carlo simulation to see how much of an effect the assumptions have.

Let us consider the general form of our game as:

If the offense always runs then the expected chance of scoring is simply \(A\).

Let us make the following assumptions:

  1. Running against a pass defense is always better than running against a run defense: \(A<B\).
  2. Passing against a run defense is always better than passing against a pass defense: \(C>D\).
  3. Passing against a run defense is always better than running against a run defense: \(A<C\).
  4. Running against a pass defense is always better than passing against a pass defense: \(B>D\).

This simple set of assumptions ensures that no strategies are dominated for either player and so our Nash equilibria will be a mixed Nash equilibrium. Let us find the expected chance of scoring for the offence at the Nash equilibrium.

The defense will ensure that the offense is indifferent:

which is equivalent to:

This has solution:

(which thanks to our assumptions is indeed a probability distribution.)

Now thanks to this we can compute \(\alpha\) which will be some measure of how much better things are for the offense at the Nash equilibria (calculated as \(1.24\) earlier on).

which we can compute at Nash equilibrium as:

Monte carlo simulation

Now that that hard work is done let us write some very simple Python code that will randomly sample values for \(A,B,C,D\) according to our assumptions and then we can see what effect this has on \(alpha\).

Let us sample our input parameters according to the following:

  • \(A \sim \mathcal{N} (60,10^2)\): normal distribution with mean 60 and standard deviation 10. Running against a run defense would work on average 60% of the time with a fair bit of variation.
  • \(B \sim \mathcal{N} (95,5^2)\): normal distribution with mean 95 and standard deviation 5. Running against a pass defense works predictably well.
  • \(C \sim \mathcal{N} (85,15^2)\): normal distribution with mean 85 and standard deviation 15. Passing against a run defense works well but not robustly.

  • \(D \sim \mathcal{N} (50,20^2)\): normal distribution with mean 50 and standard deviation 20. Passing against a pass defense does not well very well.

After sampling each parameter, a set of parameters is only accepted as valid if it follows the assumptions mentioned earlier (all the inequalities). Here is what the parameters look like after running 100000 simulations:

Now to look at some of the results. Here is what \(y\) looks like:

Recall that \(y\) is the probability with which the defense should defend the run. Finally, let us take a look at what \(\alpha\) looks like:

Recall that \(\alpha\) denotes the ratio of the scoring probability when the Nash equilibrium is played over the ‘just run all the time’. The graph is rather unhelpful as the tail is extremely long but there are some instances of our Monte Carlo simulation that have yielded a ten fold increase in scoring probability (this will occur when \(A\) is pretty low to begin with).

The mean ratio (based on our assumptions) is \(1.24\) and importantly the minimum is greater than \(1\).

All of this shows that even with a very broad relaxation of our assumptions it makes sense to at times run the ball. So game theory does in effect say that this was not ‘the worst play in super bowl history’ as it makes sense to at times pass in that situation.

TLDR: (Based on the assumptions) the coach who truly will at time pass in that situation will win the superbowl 24% more of the time than the coach who would only ever run.

Nonetheless this is ignoring a lot of the subtleties of the game itself, non more so than the fact the Patriots in fact came out in a run defense formation and the game was simply won by a great defensive play.

The code for all this can be found here.

February 15, 2015 12:00 AM

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 13, 2015

Vince Knight

Rock paper scissors lizard spock tournament

At the beginning of the week we played had one of my favourite class activity for my game theory class: the rock paper scissors lizard tournament :) In this post I will briefly go over some of the results.

If you are not familiar with Rock Paper Scissors Lizard Spock this is a good video that explains it:

This is the second time I play this in class and you can read about the first time on my old blog.

Here is how the game went (thanks to Geraint for noting everything down and Saniya for grabbing the pictures!):

(You can find the tikz code to draw that here.)

Here is a plot of the strategies played during the 1st, 2nd and 3rd rounds:

The overall strategy profile played is here:

Just like last year we are not exactly at Nash equilibria. In fact it seems that Scissors and Rock are being played a bit more often, so someone entering in to this game should respond by playing Spock (he vaporises Rock and smashes Scissors).

Here are the strategies that at some point won a duel:

and here are the losing strategies:

(All the code used to draw and analyse those plots is here.)

Hopefully my class found this interesting and fun.

February 13, 2015 12:00 AM

February 07, 2015

Vince Knight

Playing against a mixed strategy in class

This post if very late but I have been very busy with some really exciting things. I will describe some of the data gathered in class when we played against a random strategy played by a computer. I have recently helped organise a conference in Namibia teaching Python. It was a tremendous experience and I will write a post about that soon.

We played against a modified version of matching pennies which we have used quite a few times:

I wrote a sage interact that allows for a quick visualisation of a random sample from a mixed strategy. You can find the code for that at the blog post I wrote last year. You can also find a python script with all the data from this year here.

I handed out sheets of papers on which students would input their preferred strategies (‘H’ or ‘T’) whilst I sampled randomly from 3 different mixed strategies:

  1. \(\sigma_1 = (.2, .8)\)
  2. \(\sigma_1 = (.9, .1)\)
  3. \(\sigma_1 = (1/3, 2/3\)

Based on the class notation that implies that the computer was the row player and the students the column player. The sampled strategies were (we played 6 rounds for each mixed strategy):


Round 1

This mixed strategy (recall \(\sigma_1=(.2,.8)\)) implies that the computer will be mainly playing T (the second strategy equivalent to the second row), and so based on the bi-matrix it is in the students interest to play H. Here is a plot of the mixed strategy played by all the students:

The mixed strategy played was \(\sigma_2=(.71,.29)\). Note that in fact in this particular instance that actual best response is to play \(\sigma_2=(1,0)\). This will indeed maximise the expected value of:

Indeed: the above is an increasing linear function in \(x\) so the highest value is obtained when \(x=1\).

The mean score for this round by everyone was: 2.09. The theoretical mean score (when playing the best response for six consecutive games is): \(6(-.2\times 2+.8)=2.4\), so everybody was not too far off.

Here is a distribution of the scores:

We see that there were very few losers in this round, however no students obtained the best possible score: 7.

Round 2

Here the mixed strategy is \(\sigma_1=(.9,.1)\), implying that students should play T more often than H. Here is a plot of the mixed strategy played by all the students:

The mixed strategy played was \(\sigma_2=(.05,.95)\). Similarly to before this is close to the actual best response which is \((0,1)\) (due to the expected utility now being a decreasing linear function in \(x\).

The mean score for this round by everyone was: 10.69 The theoretical mean score (when playing the best response for six consecutive games is): \(6(.9\times 2-.1)=10.2\), which is less than the score obtained by the class (mainly because the random sampler did not actually pick T at any point).

Here is a distribution of the scores:

No one lost on this round and a fair few maxed out at 12.

Round 3

Here is where things get interesting. The mixed strategy played by the computer is here \(\sigma_1=(1/3,2/3)\), it is not now obvious which strategy is worth going for!

Here is the distribution played:

The mixed strategy is \(\sigma_2=(.61,.39)\) and the mean score was -.3. Here is what the distribution looked like:

It looks like we have a few more losers than winners but not by much. In fact I would suggest (because I know the theory covered in Chapter 5 of my class) that the students were in fact indifferent against this \(\sigma_1\). Indeed:


In fact, this particular \(\sigma_1\) ensures that the expected result for the students is not influenced by what they actually do:

What strategy could the students have played to ensure the same situation for the computer’s strategy? At the moment, the mixed strategy \(\sigma_2=(.61,.39)\) has expected utility for player 1 (the computer):

As this is an increasing function in \(x\) we see that the computer should in fact change \(\sigma_1\) to be \((1,0)\).

Thus if the original \(\sigma_1\) of this round is being played, so that the choice of \(\sigma_2\) does in fact not have an effect, students might as well play a strategy that ensures that the computer has no incentive to deviate (ie we are at a Nash equilibrium).

This can be calculated by solved the following linear equation:

which corresponds to:

which gives a strategy at which the computer has no incentive to deviated: \((1/2,1/2)\). Thus at \(\sigma_1=(1/3,2/3)\) and \(\sigma_2=(1/2,1/2)\) no players have an incentive to move: we are at a Nash equilibrium. This actually brings us back to another post I’ve written this term so please do go take a look at this post which involved students playing against each other and comparing to the Nash equilibrium.

February 07, 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 29, 2015

Vince Knight

Recreating golden balls in class

Here is a blog post that mirrors this post from last year. I will describe the attempt at playing Golden Balls in class.

The purpose of this was to discuss the concept of Dominance. If you are not familiar with Golden Balls take a look at this video:

This is what we did.

Firstly we played a series of four games where the score of the row player would correspond to a total of chocolates that both players would share at the end of the game (given that all the games are zero sum this is equivalent to the opposite of the score of the column player).

  1. No strategy is dominated:

    In this instance both players went for the column strategy. There is no real explanation for this: in essence they got lucky :) Thus at this stage we had 1 bar of chocolate.

  2. The first row strategy is dominated:

    From here it is obvious that the row player would go for their second strategy. This indeed happened and the column player went for their first strategy (which in fact makes no difference: the column strategy could have played either strategy). Thus at this stage we had 2 bars of chocolate.

  3. No strategies are dominated:

    This is very similar to the first game except I upped the number of chocolate bars available. Both players played their second strategy and thus we had a total of 4 bars of chocolate at this stage.

  4. The first row strategy is dominated:

    The row played went with their first strategy (as expected given the domination) however the column player went with their second strategy. This is possibly due to the slightly ‘fake’ setup of the game in terms of chocolates. Picking the second strategy ensured that they would not lose chocolates.

At the end of this I added another chocolate bar to the stash so that we had a nice even number of 6. At this point the players actually played a version of Golden Balls:

The utility in that game corresponded to the share of the chocolates: so a score of 1 implies they would both get 6, a score of 1.5 would imply they both got 9.

Both player here managed to stay away from the Nash equilibrium (which is the pair of second strategies due to iterated elimination of dominated strategies) and ended up with 6 chocolates each. Well done to G and K who were good sports and without whom we would not have been able to play the game.

This was in stark contrast with the cool result from last year.

After this I proceeded to play a best out of three game of tic-tac-toe with J to get to the idea of a game like that being pre defined: no one should really ever win that unless someone makes a mistake (or indeed plays irrationally: on that note I owe J a chocolate bar).

Here is Randall Munroe’s solution of tic-tac toe.

This leads us to the idea of best responses which is the subject of another game we played in class and one for which I’m about to start reading in the data. If you’re inpatient take a look at the corresponding post from last year.

January 29, 2015 12:00 AM

January 26, 2015

Vince Knight

Introducing Game Theory to my class

Here is a blog post that mirrors this post from last year.

I will be using my blog to extend the class meetings my Game Theory class and I have.

After playing the 2/3rds of the average game (you can see a plot of the results from last year and this year in the comments of the intro chapter).

After this we moved on to normal form games, and in particular discussed the matching pennies game:

“Two players each show a coin with either ‘Heads’ or ‘Tails’ showing. If both coins match then the 1st player (the row player) wins, otherwise the 2nd player (the column player) wins.”

This can be described as:

If you’d like to read a description of what each number represents take a look at my blog post from last year.

I asked all the students to get in to pairs and play against each other, collecting all the results with forms that you can find at this github repo.

After this we changed the game slightly to look like this:

The main point of this is to make sure that everyone understands the normal form game convention (by breaking the symmetry) and also to make it slightly more interesting (the row player now has more to win/lose by playing Heads).

I played both games with S and recorded the results to demo what was happening:

This is not what this post is about: I’m going to analyse the data collected :)

Here’s a python script that contains the data as well as the matplotlib code to obtain the plot for the initial game.

The plot shows the probability with which players played ‘Heads’ as the rounds of the game played out:

We see that we are very close at an overall probability of playing ‘Heads’ with probability 0.5. This is more or less what is expected. Now for something a bit more interesting.

Here’s a python script that contains the data as well as the matplotlib code to obtain the plot for the modified game.

This is (just like last year) not quite what we expect (which is cool). It’s pretty late as I’m writing this and I need to head to sleep so I’ll point you towards the post from last year (here) and encourage you to read that and perhaps offer your own interpretation of what is happening :)

Finally: a big thanks to my students for engaging so much, I really appreciate it.

January 26, 2015 12:00 AM

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)$ such that

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, in such a way that the above equation holds.

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

Subgroup Explorer

I’ve written a subgroup lattice generator for all groups of size up to 32. 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.

Click Go! below to refresh the viewer, or if it doesn’t load.

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 of subgroups 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.

In the image at the top of the post, the group C15 : C4 (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.

And here’s the code for a version that you can run on SageMathCloud. It allows you to input much larger groups. This was used to produce the image at the top of the post. Don’t try running it here, however, since the SageCellServer doesn’t have the database_gap package installed.

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

January 17, 2015 12:00 AM

January 16, 2015

Vince Knight

On a paper about conscientiousness

After hearing about it on TWIS I spent some time reading Other-rated personality and academic performance: Evidence and implications by Poropat. This paper is a meta analysis of various works and (TLDR): indicates that ‘intelligence’ is not as good an indicator of academic performance as is ‘conscientiousness’ (my loose interpretation of this is: ‘willingness to work hard’) and ‘openness’ (my loose interpretation of this is: ‘curiosity and interest in a subject’).

My ears really perked up when this was mentioned on the TWIS podcast as it is something I have always believed myself. This is possibly something to do with my interaction with research students. I also believe it is linked to my own educational trajectory that really benefited from my high school physics teacher who managed to show me that hard work paid off. I blogged about that here where this photo (showing a report card claiming that I did not work hard enough) is posted:

The paper

Overview of factors that potentially have an impact on academic performance.

The paper starts off by giving an overview of the ‘general intelligence factor’ denoted by g which has apparently been closely associated to academic performance. The contrast to this (again discussed in the paper) is the Five Factor Model for personality which maps individual personality to the following 5 dimensions:

  1. Agreeableness
  2. Conscientiousness
  3. Emotional Stability
  4. Extraversion
  5. Openness

One of the difficulties I had with understanding this paper was with the psychological vocabulary that I was not familiar with. Those 5 dimensions are not always easy to define but a loose (and almost certainly) incorrect interpretation of Conscientiousness is an individual’s ‘work-hard-ability’. My interpretation of Openness is an individual’s ‘want-to-learn-stuff-ability’ but the paper goes in to a pretty good discussion of each of those so I would recommend taking a look.

The paper gives a nice description and review of each dimension. I am mainly going to concentrate on what the paper says about Openness and Conscientiousness but there was one particular thing said about Extraversion that I thought was worth noting:

“Although more extraverted students may get greater attention leading to higher performance at primary level, the reduced strength of teacher student relationships at higher levels of education appear to eliminate this effect.”

Whilst this effect no longer being present at higher levels of education is perhaps a positive I feel that it does not show a ‘levelling of the playing field’ but rather that the student specific education has a lesser emphasis at higher levels of education (that is a beast of a long and terribly written sentence, it is late and I am too tired to fix it so instead added this parenthesis that makes it even longer; if you have read this far: I congratulate you). I will not dwell on this as I am not sure it is the main point of the paper nor that there is a quick solution (I teach 150+ students in my first year class…) but I thought it was interesting.

Self vs non-self evaluation

One of the important things when trying to correlate academic performance and personality is obviously getting the correct measurement for the Five Factor Model. An in depth overview of self versus non self evaluation is then given in the paper and Poropat describes how various studies have shown that non-self evaluation is a better predictor of academic performance than self evaluation (note that at this point no comparison is given to the general intelligence factor - that comes later). I think this kind of points to the idea that ‘teachers and peers know you better than you know yourself’ (or at least in terms of the Five Factor Model). It is particularly relevant to the work of Poropat’s paper as he then collects studies that look at the correlation of the Five Factor Model with academic performance: in particular only non-self evaluation studies are considered.

Using the Five Factor Model

One interesting aspect of the paper is that it emphasises that certain pedagogic approaches would be better suited to certain personalities:

“For example, discovery learning approaches help students who are higher in openness to learn while students lower on openness are aided by programmed instruction…”

I need to think about this in relation to the fact that there is other research that shows that students achieve better academic performance in an active learning environment (such as discovery learning which is apparently another term for inquiry based learning).

Hard work and curiosity are better predictors of academic performance than intellect

This is one of the nicer takeaways of the paper, by analysing 16 reports of studies that linked the five factor model to academic performance; Poropat shows that the correlation of Conscientiousness is stronger than the correlation of Openness. This is in turn stronger than previously reported correlations of the general intelligence factor.

I am sure that there will be studies and findings that report different things but I know that I’ll be using this meta analysis as a basis for further pedagogic work (in particular this paper will be very helpful for my current undergraduate research student and I: we are looking at student personality in a flipped classroom).

I know that I have always had a major preference to work with students that are (according to my personal evaluation) high on the Conscientiousness scale.

TLDR: Summary

I realise that I really have rambled in the above (the notes in my notebook are far messier still) so here are 3 bullet points I would take away from this paper:

  • There exists a 5 dimensional measure of personality;
  • Non-self evaluated versions of this measure are more accurate with regards to academic performance;
  • There is evidence that hard work is a better predictor than intellect for academic performance.

If there was one point I wish all students would believe it would be that last one. I often feel that some students believe that they have to settle for a certain level of achievement (‘because that is just how smart they are’) and this is something I have never personally been satisfied with. This is mentioned in particular in Poropat’s paper as there is evidence that personality can be modified more easily than general intelligence factor (even in older students).

It also seems evident as students of higher intellect would perhaps have been used to getting by without effort. Once things became ‘hard’ then perhaps those students who are used to working hard could indeed achieve success…

 My personal experience

To finish off this blog I thought I would throw in an xkcd style graph of my own personal ‘academic journey’ (to fully understand it, the blog post about my Physics teacher I linked to earlier is probably of interest):

(Matplotlib code here if it is of interest)

As a general summary I can certainly say that I firmly believe that any/all of my academic achievements have had very little to do with my ‘intellectual ability’ as opposed to my work ethic. Here’s a quote of Larry Bird that I really fell in love with the first time I heard it:

I’ve got a theory that if you give 100% all of the time, somehow things will work out in the end.

January 16, 2015 12:00 AM

January 07, 2015

Vince Knight

My Otis King Calculator

I was given a very cool present this Christmas period by my in laws (Check out Rachel’s fitness blog here and also this pretty incredible photo of Bryn in Antartica here). It’s an old school Otis King calulator which based on wikipedia was made sometime between 1922 and 1972.

Here it is in it’s box:

The box contains some (pretty old looking) instructions and the calculator itself:

This calculator allows you to carry out multiplication and division pretty easily thanks to the use of a logarithmic scale.

This video is a pretty good demonstration of how to work it:

If you don’t have time to look at that short video here’s the basics of trying to multiply 14 by 6. The calculator is made up of two cylinders seperated by a slider with two markers:

First we put the bottom marker on 14:

Once we have done that we align the top marker with a ‘ONE’ (there are three of them) at the top or bottom (of the upper cylinder). We do this by keeping the bottom marker aligned on the 14:

The final step is to move the top marker to 6 and read the result of the product on the bottom marker:

As can be shown in the photo the result is 84.

As you can see in the video (and or might be able to work out) it’s pretty simple to invert this process to carry out division also. For example the following photo shows the last step of 84/5= 16.8.

This is such a neat little calculator and a very cool gift. Not quite sure I’ll be letting go of my smart phone just yet but this will go really nicely on the desk in my office.

January 07, 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$

Coloring edges is almost as simple as coloring vertices, so we’ll start with that.

Generating small groups

As we’ve done in previous posts, let’s start by choosing 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’ll label by cardinality, and 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. 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 things 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

December 23, 2014

Vince Knight

Using Python and Selenium to write functional tests for a Jekyll site

Over the past six months or so I’ve become a huge Jekyll fan. In this post I’ll briefly show how to write functional tests using Selenium for a Jekyll site.

What are functional tests?

This is all extremely well described in Test Driven Development with Python. Functional tests are one aspect of test driven development (TDD). They concentrate on testing that software works as it is expected to when used by the user. TDD is the (awesome) framework in which the first thing one should do when writing code is to write a test, then check that it fails and then write code that makes it pass. Or to put it simple “Obey the Testing Goat! Do Nothing Until You Have a Test” (that is a direct quote from the book I mentioned above).

Testing takes various forms, two of which are (the hyper links there go to the corresponding Python library):

  • Doctests: this involves writing tests directly in the documentation of the code you are writing.
  • Unittests: this involves writing more robust tests in a separate script.

Let us fire up a Jekyll site (skip this if you are a jekyll regular)

This is very easy to do, after installing jekyll (see jekyll installation instructions), the following will create a base site:

$ jekyll new site_for_tests
New jekyll site installed in /Users/vince/site_for_tests.

Simply navigate to that new folder and run the following jekyll command to fire up the base site:

$ jekyll serve
Configuration file: /Users/vince/site_for_tests/_config.yml
            Source: /Users/vince/site_for_tests
       Destination: /Users/vince/site_for_tests/_site
 Auto-regeneration: disabled. Use --watch to enable.
Configuration file: /Users/vince/site_for_tests/_config.yml
    Server address:
  Server running... press ctrl-c to stop.

Then open up a browser and throw in to the url, the base site will come up:

Now we are ready to write a Selenium testing framework

Writing an initial Selenium test

So now we are going to write a test that indeed checks that the website acts like we expect. First of all let us install the Python selenium library:

$ pip install selenium

Now let us write some tests! Open up a file called and fill it with this:

Run the tests:

$ python

This will open up Firefox and (assuming all is well) return nothing in the shell.

So now let us modify the test file:

Note that I am getting ready to change the base jekyll template and build my site about writing tests. Run the tests:

$ python
File "", line 6, in <module>
    assert 'How to write tests' in browser.title  # This checks that the required title is in browser title

This time around we get an assertion error :)

Now we can go about changing our site (we are doing some TDD right now). Here is the new config file (note I have only changed the title field):

Now when we run the tests we get no assertion error:

$ python

This frees us up to write another test and then write another feature etc…

Taking things further

The above is an extremely simple example of what Selenium can do and also of how to write tests.

  1. If you know how to write unit tests but are not sure about Selenium take a look at the Selenium site (you can click on a button for Python or indeed whatever interface you would like to use). That site has a good collection of what Selenium can do (check what happens when clicking on links, checking content etc…). This is also helpful:
  2. If you are happy with Selenium but not unit tests then there are a variety of great tutorials around but to be honest I cannot recommend the Test Driven Development with Python Book enough. Harry Percival did a great job.

Here are some tests I wrote today for the site my students have put together for Code Club:

In there you can see examples of all of the above (clicking on links, checking content, checking things against a database etc…) but also the way I document the code (using what is called a ‘User Story’ which explains what a user should/would see). You can also see the way to properly ‘tear down’ the tests (so that Firefox closes).

I hope this is helpful for some: in essence you can use Selenium via Python for any site, to use it with jekyll all you need to do is have the local server running.

December 23, 2014 12:00 AM

December 21, 2014

David Horgan (quantumtetrahedron)


This work is based on the paper “Exact Computation and Asymptotic Approximations of 6j Symbols: Illustration of Their Semiclassical Limits by Mirco Ragni et al which I’ll be reviewing in my next post.

The 6j symbols tend asymptotically to Wigner dlnm functions when some angular momenta are large where θ assumes certain discrete values.







These formulas are illustrated below:




This can be modelled using sagemath.


The routine gives some great results:

For N=320, M=320, n=0, m=0, l=20, L=0,  Lmax=640

Wigner 6j vs cosθL


For N=320, M=320, n=0, m=0, l=10, L=0,  Lmax=640

Wigner 6j vs cosθL


For N=320, M=320, n=0, m=0, l=5, L=0,  Lmax=640

Wigner 6j vs cosθL



by David Horgan at December 21, 2014 11:20 PM

December 17, 2014

Martin Albrecht


When a year ends people make lists. I can only guess that several people are currently busy with writing “The 5 most revised papers on eprint ” and “The 8 best IACR flagship conference rump session presentations of 2014”. Since all the good lists are taken, my list has to be a little bit more personal. Alas, here is my list of stuff that happened in open-source computational mathematics in 2014 around me. That is, below I list what developments happened in 2014 and try to provide an outlook for 2015 (so that I can come back in a year to notice that nothing played out as planned).

If you are interested in any of the projects below feel invited to get involved. Also, if you are student and you are interested in working on one of the (bigger) projects listed below over the summer, get in touch: we could try to turn it into a Google Summer of Code 2015 project.

fplll: lattice reduction

This year I mainly worked on lattice-based cryptography. At the heart of this line of research is the assumption that finding short vectors in discrete subgroups of \mathbb{R}^n (think of a vector space, but only integer linear combinations are allowed) is hard. The main tool for finding such short vectors is lattice reduction.

The first algorithm for lattice reduction was the LLL algorithm which proved to be incredibly useful for many applications in cryptography and beyond. LLL is implemented in NTL, Pari, fplll and – as of this summer – in FLINT 2; fplll typically seems to be the fastest implementation (cf. the link above for a comparison with FLINT 2).

While LLL runs in polynomial time (read: fast) it only gives a “short vector” which is about 2^{n/4} longer than the actually shortest vector in the lattice (read: not so good). This is good enough for many applications but not good enough to, e.g. solve LWE.

To produce shorter vectors we typically employ the BKZ algorithm which is parameterised by a block size. The larger the block size the better the output but the longer it takes. BKZ is implemented in NTL and fplll, typically fplll is faster.

At AsiaCrypt 2011 Chen and Nguyen presented their results of combining various known techniques for speeding up the BKZ algorithm. They call the result BKZ 2.0. They also applied BKZ 2.0 to various benchmark problems and claimed substantial improvements over the public state of the art. However, for some reasons only known to the authors and the AsiaCrypt 2011 programme committee, their paper was accepted without them publishing their source code. Since then we’re in the somewhat strange situation that everybody believes BKZ can be made to run much faster (you might even get your paper rejected because you are not using the state of the art, i.e. BKZ 2.0) but no one has publicly reproduced these results.

In autumn we took some first steps towards fixing this situation by adding better pre-processing and an easier linear pruning interface to fplll. I know that others have patched fplll as well, but as far as I know they never made their changes public. In the process, we also moved fplll’s development to Github, so send your pull requests.

While the result is an improvement over plain BKZ, there’s still a lot of work to be done to come even close to the results of Chen and Nguyen. For example, they use extreme pruning instead of the simple linear pruning strategy currently implemented in fplll and it’s not clear how to pick pruning parameters for extreme pruning … however, there’s a paper on the arXiv which promises an answer to this question. Furthermore, currently the user has to pick pre-processing parameters by hand, something the implementation should take care off by default. Finally, a recent paper on ePrint claims that phase-based enumeration is much faster than the Kannan-Helfrich enumeration algorithm which is implemented everywhere including fplll. I’d consider addressing any of these issues a valuable contribution.

dgs: sampling from a discrete Gaussian distribution

A central step in most lattice-based cryptography is to sample from a discrete Gaussian. A discrete Gaussian distribution over the Integers is a distribution where the integer x is sampled with probability proportional to \mbox{exp}(-(x-c)^2/(2\sigma^2)) where σ is the width parameter (close to the standard deviation) and c is the centre. There are several algorithms to choose from to sample from such a distribution each being better or worse in some situations. Somewhat surprisingly, though, when I got started implementing some lattice-based crypto there was no C library available that allows to sample from a discrete Gaussian with reasonable efficiency. Now there is. The library is called dgs and it is included in Sage by default and also unpins our GGHLite implementation, so it has seen some usage.

A few things still need to be done, though. Some bugs were fixed in the stand-alone library but not ported back to Sage, yet. Also, the library would benefit from more tests being run by make check. Finally, implementing the Discrete Ziggurat algorithm would complete the picture.

oz: computing in some Cyclotomic number rings

Most of the interesting code in our GGHLite implementation is in the oz submodule which implements arithmetic in \mathbb{Z}[x]/(x^n+1), where n is a power of two, and with all operations in quasi-linear time. While the code is already fairly modular, i.e. we separated crypto applications from arithmetic, it might make sense to outsource this module into a separate library so it can used more easily by other projects should they wish to.

When we are doing this, we should probably also split up the dgsl module. This module implements sampling from a discrete Gaussian distribution over arbitrary lattices (in contrast to dgs which implements it over the integers). This module contains two essentially independent parts. One part samples from lattices represented by a basis matrix using the GPV algorithm. Another part samples from ideal lattices represented by an ideal generator in \mathbb{Z}[x]/(x^n+1) using Peikert’s algorithm. The latter relies heavily on oz (as well as dgs) and might as well be moved there, the former has no connection to oz and could be either included in dgs (which would entail making FLINT a dependency of dgs) or remain independent.

Linear algebra over small finite fields

I didn’t work as much on linear algebra over small finite fields as I would have liked to in 2014. I doubt I’ll make it a priority of mine in 2015 either, so if anybody wants to jump in to help, that’d be much appreciated.


M4RI implements dense linear algebra over \mathbb{F}_2, i.e. the field with two elements 0 and 1. We released one bugfix release of M4RI this year.

If you read this blog, you probably know that Enrico’s implementation of Gaussian elimination is faster than our own. As far as I can tell the advantage of GF2Toolkit over M4RI comes from avoiding a lot of management overhead. To illustrate this point, consider the following output of Google’s perf tools on running M4RI’s Gaussian elimination on a 4096 x 4096 dense full rank matrix:

     .      . 1121:      switch(__M4RI_M4RM_NTABLES) {
    42    131 1122:   case 8: t[7] = T[ 7]->rows[ L[7][ (a >> 7*k) & bm ] ];
    57     79 1123:   case 7: t[6] = T[ 6]->rows[ L[6][ (a >> 6*k) & bm ] ];
    27     47 1124:   case 6: t[5] = T[ 5]->rows[ L[5][ (a >> 5*k) & bm ] ];
    14     25 1125:   case 5: t[4] = T[ 4]->rows[ L[4][ (a >> 4*k) & bm ] ];
    29     52 1126:   case 4: t[3] = T[ 3]->rows[ L[3][ (a >> 3*k) & bm ] ];
    15     34 1127:   case 3: t[2] = T[ 2]->rows[ L[2][ (a >> 2*k) & bm ] ];
    14     27 1128:   case 2: t[1] = T[ 1]->rows[ L[1][ (a >> 1*k) & bm ] ];
     6     17 1129:   case 1: t[0] = T[ 0]->rows[ L[0][ (a >> 0*k) & bm ] ];
     .      . 1130:         break;
     .      . 1131:   default:

     .      . 1137:   switch(__M4RI_M4RM_NTABLES) {
   970   1946 1138:     case 8: _mzd_combine_8(c, t, wide); break;

The numbers in the first two columns indicate how much time we spent in each line. As you can see, we’re spending between 20% and 25% of the time it takes to perform additions (_mzd_combine_8) with setting them up (everything else); We are performing 4096/128 · 8 = 256 XORs in _mzd_combine_8, which isn’t much and so our setup overhead is hurting us.


M4RIE is a library for fast arithmetic with matrices over small even characteristic fields, i.e. \mathbb{F}_2[x]/f(x) where f(x) is an irreducible polynomial over \mathbb{F}_2[x] of degree up to 16. We released one bugfix release of M4RIE this year.

Last year I added some code to M4RIE for computing with matrices over \mathbb{F}_2[x], i.e. where the entries are high(er) degree polynomials over \mathbb{F}_2. The strategy I implemented is some “evaluation, pointwise-multiplication, interpolation” scheme where I use Dan Bernstein’s “Optimizing linear maps modulo 2” to cut down the cost of first and last step. Unfortunately, I didn’t get around to work more on this code this year. While I still don’t know an application for this, it would be fun to see how far we can push this. But I guess to do this properly, we’d need to also take another look at the Number Theoretic Transform to realise such multiplications, at least when the dimension of the matrices is not much bigger than the degree of the polynomials.

Another area for improvement is that the formulas we use to realise multiplication for degrees up to 16 are not always optional. In particular, we know that the following improvements are possible for degree 6 (18 → 15), degree 8 (27 → 24), degree 9 (31 → 30), degree 10 (36 → 33), degree 11 (40 → 39), degree 12 (45 → 42), degree 13 (49 → 38), degree 14 (55 → 51), degree 15 (60 → 54), degree 16 (64 → 60), where the numbers in brackets are the current and the best known number of multiplications over \mathbb{F}_2. Some of these improvements can be realised by simply dropping in known better formulas, some of them would be a bit more involved because they rely on finite field embeddings.


CryptoMiniSAT is the SAT solver with the best integration into Sage. However, Sage is still using CryptoMiniSAT 2.9.6 instead the more recent 4.x series. This is partly because we can’t get our act together and partly because Máté decided to go with CMake instead of Autotools in the current CryptoMiniSat series. There’s a pull request idling around which improves Autotools support for CryptoMiniSAT. I should probably follow up on this. Once this is taken care of (which shouldn’t take more than 1-2 hours), we should update CryptoMiniSAT in Sage. The interface of CryptoMiniSAT hasn’t changed much, so this second step shouldn’t be too hard either. Some options might have changed, so I would guess solverconf_helper.cpp would need to be adapted slightly.


Sage saw four major releases in 2014. Sage 6.1 in February, Sage 6.2 in May, Sage 6.3 in August and Sage 6.4 in November.

My main contributions this year were to add better support for computing with lattices: I mainly worked on the fplll interface (see above) and on sampling from discrete Gaussian distributions (also above). I also fixed the occational bug and reviewed the occational ticket, but not as much as I would have liked to. In fact, there’s just been a friendly reminder that we have way too many tickets for Sage with status “needs review” which means that someone contributed some code and that code is now waiting to be reviewed so it can be included (or revised).

As always there’s much to be done for Sage, but too little time. Here are some examples.

SAT Solvers

Besides the tasks listed under the CryptoMiniSAT heading, more work should be done on Sage’s SAT solving capabilities. Currently, Sage will fail to solve any SAT problem without installing additional software because it believes that no SAT solver is included by default. Turns out, this is not correct. That is, Sage ships with GLPK and GLPK can be used as a SAT solver. That won’t break any performance records but it’s better than nothing. So we should use GLPK as the poor person’s SAT Solver. Indeed, adding support for this should be fairly straight forward. Here’s what Nathann Cohen had to say who pointed me towards GLPK:

There are some sentences about it at the end of that :

And this page says that there is a dedicated pdf in GLPK’s doc:

But it seems that the interface is pretty basic, and may have to work through files… or though LP ! But as you can already produce DIMACS SAT instances with you code, perhaps you can just call GLPK on that? It would be better than nothing, plus you can say that the feature is standard, and also write everywhere that users can download a “real solver” if they want to.

Secondly, while Sage currently is happy to write DIMACS files (the standard format for SAT problems) it does not read DIMACS files.


For a while now I’ve had adding an interface to FGB on my TODO list. FGB is Jean-Charles Faugère’s implementation of the F4 Gröbner basis algorithm. It is not open source, but it is a good implementation of F4, something which isn’t exactly widely available in the open-source world. Besides, being able to easily compare with FGB surely would be useful.


When I worked on solving multivariate systems of equations with noise I added support for the SCIP constraint integer programming solver to Sage. Constrained integer programming allows to solve non-linear systems of equations and inequalities. For example, it allows to model and solve systems which contain constrains like x \cdot y + 4x + 1 \geq 0 and z = x\ \mbox{OR}\ y which comes in handy from time to time. SCIP is open source in the sense that you can read the source code, but its license does not live up the standards of the Open Source Initiative and we can’t ship it. Still, it is typically faster than real open-source solutions and the developers are happy to help. I’ve not touched this code in a while because I don’t work in this area at the moment, so someone who does should pick up the project.

by martinralbrecht at December 17, 2014 04:53 PM

Liang Ze

Lattice of Subgroups II - Coloring Vertices

In my previous post, I showed how to use Sage to generate the subgroup lattice of a group, and define labels for the subgroups. In this post, I’ll demonstrate how to color subgroups with different colors according to some desired property. If you’re not interested in code, scroll to the bottom of the post for a visual collection of groups and their subgroup lattices.

Lattice of the dicyclic group $Dic_3$

First, let’s rerun the code from the previous post. We’ll choose a group we like and generate its poset. For this post, I’ll label the subgroups by their cardinality. If you’re trying this code in SageMathCloud or your own version of Sage that has the database_gap package installed, I strongly recommend labelling the subgroups using structure_description() instead).

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

Normal subgroups

Now suppose we wish to know which subgroups are normal. We can do so with the following:

Coloring takes place after labelling, so when we’re defining the color dictionary, we have to use the (re)labelled vertices rather than the original vertices. Hence the use of label[x] instead of just x.

Here are some more examples. They’re all just variations of the above.

Abelian subgroups and the center

Let’s say we want to highlight the abelian subgroups, with special emphasis on the center of the group. We can do this with a slight modification to the color dictionary:

Maximal subgroups and the Frattini subgroup

The Frattini subgroup is the intersection of all the maximal subgroups of $G$. We can see this by highlighting the Frattini and maximal subgroups.

Sage has a built-in function for getting the Frattini subgroup. To get the maximal subgroups, however, we’ll have to find the elements covered by the greatest element (or top) of the poset $P$.

Sylow subgroups

Getting the Sylow $p$-subgroups takes a little more work, since Sage doesn’t have a single function that generates all the Sylow subgroups at once.

In Sage, G.sylow_subgroup(p) returns one Sylow $p$-subgroups. To get all the Sylow $p$-subgroups, we could take all conjugates of this Sylow subgroup (since all Sylow $p$-subgroups are conjugate). A faster way, however, is to use the fact that the cardinality of all Sylow $p$-subgroups is the maximal $p^{th}$ power dividing the order of $G$.

More examples

In the next post, we’ll look at labelling edges. This is particularly useful if we want to determine if $G$ has subgroup series with certain properties.

December 17, 2014 12:00 AM

December 15, 2014

Liang Ze

Lattice of Subgroups

This is the first in a series of posts on visualizing groups via their lattice of subgroups.

Lattice of the dihedral group $D_4$

Displaying the Lattice of Subgroups

One way of getting a better understanding of a group is by considering its subgroups. The lattice of subgroups (more precisely, the Hasse diagram of this lattice) gives us a way to visualize how these subgroups relate to each other and to their parent group. Here’s how to do it in Sage:

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

By default, the vertex labels of the Hasse diagram will be the description of the object that the vertex represents. In our case, something like Subgroup of (Dihedral group of order 8 as a permutation group) generated by [(1,2)(3,4)], which would be way too long to display nicely. Because of this, in the above code, I’ve decided not to label the vertices. I’ve also chosen to make the vertices hexagonal and white.

Relabelling vertices

Without labels, it can be hard to tell what subgroups we’re looking at. We can define new labels for these vertices by defining a dictionary where the keys are the original vertices and their corresponding values are our new labels.

Labelling by generators

One way to tell what the subgroups are is to look at their generators:

This isn’t very pretty, and just knowing the generators doesn’t give us much intuition about the group.

Labelling by cardinalities

Alternatively, we could label the subgroups by their cardinalities:

If you ran the preceding code, you probably encountered an error message. This is because Sage currently requires that vertex labels be injective i.e. distinct vertices must have distinct labels. There’s a quick but slightly ugly fix for this: just pad spaces around the labels to make them all unique:

Labelling by structure description

However, cardinalities still don’t tell me very much about the subgroup. Fortunately, Sage has a method for describing the structure of a small group: H.structure_description() where H is the group in question.

Unfortunately, this method requires the GAP group database, which is not installed with the Sagecell version of Sage. However, the free SageMathCloud service’s installation of Sage does have the group database installed, so you can try the following code there. This code was used to produce the image at the start of this post:

More examples

Try playing around with different groups and different labelling methods!

And here are some questions that might arise while playing around with subgroup lattices:

  • In the code, I’ve technically only defined the poset of subgroups. However, it turns out the poset of subgroups is also always a lattice. Why?
  • When is the subgroup lattice also distributive?
  • When is the subgroup lattice a chain?

In the next post, we’ll add some color to the subgroup of lattices by coloring the subgroups according to properties they have.

December 15, 2014 12:00 AM

December 14, 2014

Vince Knight

A busy term

Here’s my final reflective post for Imogen Dunne’s final year project (you can find the first here, the second here and the third here).

The main feeling I have as I reflect on this past term is how tired I am.

This term/semester/thing has been so amazingly busy. Before I go any further I need to thank my PhD student Jason Young. Jason has just started his PhD and has been acting as a full teaching assistant for me with this course. I have no idea how this course ran without having someone to help me last year and this leads me to my first point of reflection: students have engaged very well with the course.

I believe (hope is perhaps a less biased word) that the main reason this has been busy is because this has been an active learning experience for my students. This could be due to the flipped learning environment (I would like to think so) but also perhaps just having a really greatly engaged set of students.

Here is a list of things that I plan on changing for next year:

  • Another set of videos - I’m going to double the amount of videos I have. I think this is a natural thing to do as I’m more aware of the difficulties students have. This is all part of the reflective pedagogy that revolves around the concept of being able to react in a timely manner to feedback as to student difficulties (if I’m not doing this than I am no better than a book). The reflective approach ensures a dynamic reaction to difficulties which can happen on various scales:

    • Class meetings: based on difficulties during the week (micro level)
    • From year to year: based on major trends during the year (macro level)
  • More scaffolding on technical issues - I always underestimate the starting point of some students. I need to do a better job helping them with simple things like using a mouse, using internet browsers and also more tricky things like debugging LaTeX.

  • As discussed before: a much better scaffolding of student tutors (the undergrads).

On the whole though I am very pleased with how this year has gone. In particular I feel that a number of students have not only learnt to code but also understood the importance of learning to code in conjunction with learning mathematics. Here is a quote that I’m taking from one of the 3 page papers that students have had to hand in at the end of term (this particular one looking at the futurama theorem):

“It has to be said, when I first looked at Keeler’s proof, the whole thing did seem incredibly complicated. But after coding it and using it in ‘real life’, it is pretty straight forward.”

This is really nice to read as it’s something I say a lot: coding complicated mathematics will often help understand it.

It reminds me a bit about this nice lightning talk that one of the students gave last year:

December 14, 2014 12:00 AM

December 11, 2014

David Horgan (quantumtetrahedron)


In order to follow up some work on the the last two posts I have been looking at  the capabilities of Sagemath with regard to calculating hypergeometric functions:

Fortunately sagemath can implement at number of great codes for hypergeometric functions including:

These give excellent results:





The ability to use mpmath code is very useful since it enables me to calculate a wide range of hypergeometric functions as can be seen on the reference pages. I’ll be using this over the coming posts starting with :

Exact Computation and Asymptotic Approximations of 6j Symbols:
Illustration of Their Semiclassical Limits which is in preparation.

Related articles


by David Horgan at December 11, 2014 08:49 PM

December 10, 2014

Vince Knight

A Sneak Preview of Game Theory in Sage (3/3): Normal Form Games

In two previous posts I have discussed two game theoretical concepts that are in/on their way in to Sage:

Since writing those and alluding to more development that myself and an undergraduate here at Cardiff are working on, I’ve had a fair few people asking about when Normal Form Games will be included…

The purpose of this post is to say: extremely soon! A NormalFormGame class has now also been merged in to the develop branch of Sage (so it will be in the next release).

What is a normal form game?

These are sometimes referred to as bi-matrix games or strategic form games. I wrote a blog post about these in reference to choosing a side of the pavement to walk on: here.

In this post I’ll take a look at what the new Sage class allows you to do.

Consider the game I used to model two individuals walking on either the left or the right of the pavement:

\[ A = \begin{pmatrix} 1&-1\\
-1&1 \end{pmatrix} \] \[ B = \begin{pmatrix} 1&-1\\
-1&1 \end{pmatrix} \]

Matrix \(A\) gives the utility to the first person (assuming they control the rows) and the matrix \(B\) gives the utility to the second person (assuming they control the columns). So if both individuals walk on their left then they both get a utility of 1 (ie they don’t bump in to each other).

Defining a game

We can define these two matrices in Sage and will be able to define a NormalFormGame as follows:

sage: A = matrix([[1, -1],[-1, 1]])
sage: B = matrix([[1, -1],[-1, 1]])
sage: g = NormalFormGame([A, B])

As you can see the NormalFormGame class uses a list of two matrices to construct a game. If you look at the documentation you’ll see that there are other ways to construct games. To see that this has indeed worked we can just see the output of the game:

sage: g
Normal Form Game with the following utilities: {(0, 1): [-1, -1], (1, 0): [-1, -1], (0, 0): [1, 1], (1, 1): [1, 1]}

This displays a dictionary of the strategy:utility pairs. The form of the output is actually based on gambit syntax.

We can use this class to very easily obtain equilibria of games:

Finding Nash equilibria

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

The output shows that there are three Nash equilibria for this game:

  • Both players walking on their right;
  • Both players walking on their left;
  • Both players alternating from left to right with 50% probability.

There are currently 2 algorithms implemented in Sage to calculate equilibria:

  1. A support enumeration algorithm (this is not terribly efficient but is written in pure Sage so will work even if optional packages are not installed and for typical game sizes will work just fine).
  2. A reverse search algorthim which calls the optional ‘lrs’ library.

My student and I are currently actively developing further integration with gambit which will allow for a linear complementarity algorithm and also solution algorithms for games with more than 2 players.

Here is one other example:

sage: A = matrix([[0, -1, 1, 1, -1],
....:             [1, 0, -1, -1, 1],
....:             [-1, 1, 0, 1 , -1],
....:             [-1, 1, -1, 0, 1],
....:             [1, -1, 1, -1, 0]])
sage: g = NormalFormGame([A])
sage: g.obtain_nash(algorithm='lrs')
[[(1/5, 1/5, 1/5, 1/5, 1/5), (1/5, 1/5, 1/5, 1/5, 1/5)]]

As you can see above, I’ve created a game by just passing a single matrix: this automatically creates a zero sum game. I’ve also told Sage to make sure it uses the 'lrs' algorithm (although 'enumeration' would handle this 5 by 5 game just fine).

Finally if you’re not actually sure what that game is take a look at this little video:

I’m very excited to see this in Sage (soon!) and am actively working on various other things that I know at least I will find useful in my research and teaching but hopefully others will also.

December 10, 2014 12:00 AM

November 22, 2014

Vince Knight

Thinking about divisibility by 11

This post is based on a class meeting I had recently with my programming class. It was based on trying to use code to help identify a condition that a number must obey for it to be divisible by 11. Readers of this blog might be aware that the following is incorrect but stick with me.

Exploring a statement

A number is divisible by 11 if and only if the alternating (in sign) sum of the number’s digits is 0.

To help with notation let us define \(f:x\to\text{alternating sum of digits of x}\) so for example we have:


It is immediate to note that for \(N< 100\) \(f(N)=0\) if and only if 11 divides \(N\) (\(11\;|\;N\) for short). Before trying to prove our statement we could check it for a few more numbers:

  • \(f(110)=0\)
  • \(f(121)=0\)
  • \(f(132)=0\)
  • \(f(143)=0\)
  • \(f(154)=0\)
  • \(f(165)=0\)
  • \(f(176)=0\)
  • \(f(187)=0\)
  • \(f(198)=0\)

So things are looking good! We could now rush off and try to prove that our statement is correct… or we could try more numbers. The easiest way to ‘try enough’ is to write some simple code (the following is written in Python):

class Experiment():
    A class for an experiment
    def __init__(self, N):
    Initialisation method:
    inputs: N - the number for which we will check the conjecture
        self.N = N
        self.divisible_by_11 = N % 11 == 0
        self.sum_of_consecutive_digits = sum([(-1) ** d *int(str(N)[d]) for d in range(len(str(N)))])
    def test_statement(self):
        Returns True if 'A number is divisible by 11 iff the alternating sum digits is 0' holds for this particular number.
        if self.divisible_by_11:
            return self.sum_of_consecutive_digits == 0
        return self.sum_of_consecutive_digits != 0

This creates a class for an Experiment for a given number, which has a couple of attributes relevant to what we’re trying to do:

>>> N = Experiment(121)
>>> N.divisible_by_11
>>> N.sum_of_consecutive_digits

There is also a method that checks the if and only if condition of our statement:

        if self.divisible_by_11:
            return self.sum_of_consecutive_digits == 0
        return self.sum_of_consecutive_digits != 0

So if the number is divisible by 11 then the statement is true if the sum is 0. If the number is however not divisible by 11 then the statement is true if the sum is not 0.

We can thus check for any given number if our statement is true:

>>> N = Experiment(121)
>>> N.test_satement()  # 121 is divisible by 11 and 1-2+1==0
>>> N = Experiment(122)
>>> N.test_statement()  # 122 is not divisible by 11 and 1-2+2!=0

So before attempting to prove anything algebraically let’s just check that it holds for the first 10000 numbers:

>>> all(Experiment(N).test_statement() for N in range(10001))

Disaster! It looks like our statement is not quite right!

The following might help us identify where (outputting a list of numbers for which the statement is false):

>>> [N for N in range(10001) if not Experiment(N).test_statement()]
[209, 308, 319, 407, 418, 429, 506, 517, 528, 539, 605, 616, 627, 638, 649, 704, 715, 726, 737, 748, 759, 803, 814, 825, 836, 847, 858, 869, 902, 913, 924, 935, 946, 957, 968, 979, 1309, 1408, 1419, 1507, 1518, 1529, 1606, 1617, 1628, 1639, 1705, 1716, 1727, 1738, 1749, 1804, 1815, 1826, 1837, 1848, 1859, 1903, 1914, 1925, 1936, 1947, 1958, 1969, 2090, 2409, 2508, 2519, 2607, 2618, 2629, 2706, 2717, 2728, 2739, 2805, 2816, 2827, 2838, 2849, 2904, 2915, 2926, 2937, 2948, 2959, 3080, 3091, 3190, 3509, 3608, 3619, 3707, 3718, 3729, 3806, 3817, 3828, 3839, 3905, 3916, 3927, 3938, 3949, 4070, 4081, 4092, 4180, 4191, 4290, 4609, 4708, 4719, 4807, 4818, 4829, 4906, 4917, 4928, 4939, 5060, 5071, 5082, 5093, 5170, 5181, 5192, 5280, 5291, 5390, 5709, 5808, 5819, 5907, 5918, 5929, 6050, 6061, 6072, 6083, 6094, 6160, 6171, 6182, 6193, 6270, 6281, 6292, 6380, 6391, 6490, 6809, 6908, 6919, 7040, 7051, 7062, 7073, 7084, 7095, 7150, 7161, 7172, 7183, 7194, 7260, 7271, 7282, 7293, 7370, 7381, 7392, 7480, 7491, 7590, 7909, 8030, 8041, 8052, 8063, 8074, 8085, 8096, 8140, 8151, 8162, 8173, 8184, 8195, 8250, 8261, 8272, 8283, 8294, 8360, 8371, 8382, 8393, 8470, 8481, 8492, 8580, 8591, 8690, 9020, 9031, 9042, 9053, 9064, 9075, 9086, 9097, 9130, 9141, 9152, 9163, 9174, 9185, 9196, 9240, 9251, 9262, 9273, 9284, 9295, 9350, 9361, 9372, 9383, 9394, 9460, 9471, 9482, 9493, 9570, 9581, 9592, 9680, 9691, 9790]

The first of those numbers is \(209=11\times19\) so it is divisible by 11 but \(f(209)=2-0+9=11\) and if we calculate \(f\) for a few more of the numbers in the above list we again get \(11\). At this point in time it seems like we need to adjust our statement.

Sufficient evidence for a conjecture

A number is divisible by 11 if and only if the alternating (in sign) sum of the number’s digits is divisible by 11.

A slight tweak of the Experiment code above gives:

class Experiment():
    A class for an experiment
    def __init__(self, N):
    Initialisation method:
    inputs: N - the number for which we will check the conjecture
        self.N = N
        self.divisible_by_11 = N % 11 == 0
        self.sum_of_consecutive_digits = sum([(-1) ** d *int(str(N)[d]) for d in range(len(str(N)))])
    def test_conjecture(self):
        Returns True if 'A number is divisible by 11 iff the alternating sum digits is 0' holds for this particular number.
        if self.divisible_by_11:
            return self.sum_of_consecutive_digits % 11 == 0
        return self.sum_of_consecutive_digits % 11 != 0

Now let us check the first 100,000 numbers:

>>> all(Experiment(N).test_conjecture() for N in range(100001))

When we have a lot of evidence for a mathematical statement we can (generally) start calling it a conjecture. At this point we probably can attempt to prove that the conjecture is true:


Let \(n_i\) be the \(i\)th digit of the \(m\) digit number \(N\), so we have \(N=\sum_{i=1}^{m}n_i10^{i-1}\). Using arithmetic modulo \(11\) we have:



The right hand side of that is of course just \(f(N)\) so \(11\;|\;N\) if and only iff \(11\;|\;f(N)\) (as required).

This is how a lot of mathematics gets done nowadays. Statements get made, then refined then checked and then finally (hopefully) proved. A nice book that describes a conjecture that stayed a conjecture for a long time (until ultimately being proved) is Proofs and Confirmations: The Story of the Alternating-Sign Matrix Conjecture by Bressoud.

November 22, 2014 12:00 AM

November 19, 2014

Liang Ze

The Argument Principle

This post illustrates the Argument Principle.

Let $f$ be a polynomial, and $C$ be (an arc of) a circle of radius $R$ centered at the origin.

The code below generates 2 plots. The first plot shows the domain of $f$. We plot the roots of $f$ together with $C$. Let $n_1$ be the number of roots contained within $C$.

The second plot shows the range of $f$. We plot $f(C)$, along with a marker at the origin. Let $n_2$ be the number of times the curve winds around the origin.

You can verify that $n_1 = n_2$. As you vary the radius $R$, observe how $C$ and $f(C)$ change, and how this affects $n_1$ and $n_2$.

November 19, 2014 12:00 AM

November 16, 2014

Martin Albrecht


Sage 6.4 is out. Here are some (to me) particularly exciting changes.

#16479: Vincent Delecroix: package for pip the Python installer

Sage now ships with pip which means you can use pip to install your favourite Python packages into the Sage environment, making it a lot easier to get access to the Python ecosystem. For example, say you want to use BatzenCA in Sage. You’d call

$ sage -pip install batzenca

and you are done.

Update: Sorry, I mispoke: pip is an optional package at the moment, not a standard package. You’ll have to sage -i pip it first.

#15915: Martin Albrecht: add discrete Gaussian samplers to Sage

I contributed code to sample from discrete Gaussian distributions to Sage. Discrete Gaussian distributions over the Integers sample integers proportionally to exp(-(x-c)²/(2σ²)), where c is the center and σ the Gaussian width parameter (roughly, the standard deviation). Here’s an example:

sage: from sage.stats.distributions.discrete_gaussian_integer import DiscreteGaussianDistributionIntegerSampler as DiscGauss
sage: D = DiscGauss(3.0)
sage: D() # output random by definition

Discrete Gaussian distributions are commonly used in lattice-based cryptography. For example, our GGHLite implementation makes use of the same code. That is, most of the code is written in C99 and also available as the stand-alone dgs library under a 2-clause BSD library.

I also implemented discrete Gaussians over arbitrary lattices. A discrete Gaussians over a lattice Λ is a distribution where point x occurs with probability:

\mbox{exp}(-|x-c|_2^2/(2\sigma^2))/(\sum_{x \in \Lambda} \mbox{exp}(-|x|_2^2/(2\sigma^2)))

However, this code not part of the dgs library but a pure Python module. A C99 implementation of the same algorithms as well as specialised algorithms for ideal lattices is part of the GGHLite implementation as the dgsl module, though.

#16803: Marc Masdeu: Reimplement matrix_integer_dense using FLINT

Marc switched over matrices over \mathbb{Z} to using FLINT as the native data type. Previously, we had our own data type (arrays of arrays of mpz_t), custom code and lots of conversions to other data types (LinBox, IML, …) to realise some functionality such as characteristic polynomials. FLINT is a sane default data type as it is quite fast for many operations. Conversions to LinBox et al. are still in place so functionality should not be lost. However, some decisions about defaults might be outdated now so report any issues on trac.

#16996: Volker Braun: IPython notebook with Sage Extensions

Thanks to Volker’s efforts you can now use the iPython notebook for Sage. To me the iPython notebook feels a lot more modern than the Sage notebook, partly because it is actively being developed whereas not much work is currently done on the native Sage notebook (it has seen some work done on it in the current release, though). I also like that iPython notebooks are files in your current directly which means you can keep them project local (others disagree about this design decision). To start the iPython notebook run:

$ sage -notebook=ipython

by martinralbrecht at November 16, 2014 08:46 PM

November 15, 2014

Liang Ze

Partitions and Posets

In this post, we generate the Hasse diagram of a set partition.

Partitions of a 4 element set

The following piece of code generates the image above.

Read on for an explanation of the code.


A partition of a set $X$ is a collection $p$ of non-empty subsets of $X$ such that $X$ is the disjoint union of these sets.

In SAGE, you can get the set of all partitions of the $N$ element set {$1,2,\dots,N$} using SetPartitions:


Each item in the preceding list is a partition of $X$. The elements of each partition are called blocks. At the top of the list, we see the trivial partition consisting of just one block, $X$. At the other end, we see the singleton partition consisting of $|X|$ blocks, where each block contains a single element of $X$. All other partitions of $X$ fall somewhere in-between these two partitions.

We can make this notion of “in-betweenness” precise by defining a relation on the partitions of $X$. We say that $q$ is a refinement of $p$ if each block of $q$ is contained in some block of $p$.

In Sage, you can see the refinements of a partition using the method refinements():


For a fixed $X$, the set $\mathcal{P}$ of all partitions of $X$ has the structure of a poset given by $q \leq p$ if $q$ is a refinement of $p$.

In Sage, we can construct a Poset by specifying an underlying set $P$ along with a function $f:P\times P \to$ {$\text{True},\text{False}$} where

The resulting poset can be visualized via its Hasse diagram, which is a directed graph with paths from $q \to p$ if $q \leq p$. We can generate a Hasse diagram of a poset using the show() method.

The final piece of code (at the top of the page) combines everything above to produce the Hasse diagram of a set. The function Partition_Poset first generates the set of partitions of an $N$ element set, then converts it to a poset. The function p_label relabels the partitions so that they look prettier. I’ve also tweaked some options in the show() method to make things look nicer.

November 15, 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 ( at November 14, 2014 02:31 PM

Vince Knight

Scaffolding tutors and how to better prepare for different pedagogies

Here’s my third reflective post for Imogen Dunne’s final year project (you can find the first here and the second here.

The course is now on the final straight as students have finished the class test which is always a bit of a milestone as they now begin to really individualise their learning experience (working on individually chosen projects etc…). As we get to this point I’ve got 3 specific things I’ve been thinking about:

  1. Class test performance

    With over half of the scripts marked it seems like students have done quite well on what is a difficult piece of assessment. It seems that the performance is overall better than last year but it’s too early to try and guess as to why that is.

  2. Scaffolding undergraduate tutors

    I use second year students to tutor in the lab sessions. They are all students who did the course last year. I don’t think I’ve done the best job explaining the roles of the tutors and this is something I need to get right next year. I still am very excited about using undergraduate tutors as I think it’s a brilliant experience for them as it continues their learning. The other (super important reason) is that I think this peer level of instruction is undoubtedly of benefit to the students on the course (if this wasn’t a quick rough drop of thoughts I’d be able to find a large quantity of resources relevant to this).

    The flipped approach used requires the tutoring to be very light touch, and more about giving feedback than giving ‘help’. I didn’t do the best job of explaining this as some undergraduate tutors felt it was their job to ‘get students code to work’. I’ll probably run a bit more of an ‘expansive’ training session with closer shadowing at the start.

    All the tutors have done a brilliant job and I’m very pleased, I just have to think carefully about how to best ‘scaffold’ and support them.

  3. Communicating the different pedagogy

    Some students were questioning the timing of my office hours: they are after the class meeting (which is in turn after the lab sessions). Whilst most students are on board and understand I think that the fact that a couple of students didn’t understand this shows that I didn’t do the best job of communicating the ideas.

    Involving the students in a discussion about the pedagogy they are experiencing is something I think is very important (especially when using something they might not have experienced before). As such I have tried and continue to try to talk about this throughout my interactions with students (from 1 to 1 meetings all the way to class meetings) but perhaps I could do more. Perhaps even just showing a picture like this would be helpful:

    That shows the reason why I have my class test at the end of a week, the idea being that most students will be happy with content through labs, some will be happy after the class meeting and a small amount will require specific time with me to help them through. This goes back to the premise of a flipped environment which aims to make best use of contact time: my class meetings are meant for us to further understanding and go towards the tip of Bloom’s taxonomy but also address specific difficulties.

On a whole I’m very happy with how this class is going, I feel that most students are engaging fully with the course and seem to be enjoying it. I also feel that I’ve gotten certain things a bit better than last year this year which could be an explanation of the slightly better mark in the class test (I still have 70 scripts to mark this weekend so it’s probably too early to tell).

One particularly nice piece of feedback is how students like the new class website (jekyll for the win!) and in particular have enjoyed the use of a comment system on all pages: it’s always nice to have that permanent record of discussions I’m having with the students which in turn could help other students. Ideally the discussion would be peer to peer but that hasn’t happened much this year (mainly me answering queries) but it has happened once which I’m happy about (small steps, big wins) :)

November 14, 2014 12:00 AM

Liang Ze


This blog is powered by the Sage Cell Server. You can type Sage/Python code into the cell below, and press Shift+Enter to evaluate it (or click “Evaluate”).

November 14, 2014 12:00 AM

October 25, 2014

Vince Knight

Busy office hours

Here’s my second reflective post for Imogen Dunne’s final year project (you can find the first here).

We are now at the end of Week 4 of the course and I’m glad to say that I think overall things are going well.

  • My office hours have gotten very busy. This is awesome. Students are coming to see me, genuinely having struggled with concepts and this is often a result of myself or another tutor identifying specific issues in a lab session and saying something like:

    I’ll give you the tick for that but come and see me during office hours to talk about it as I think you’re a bit confused there.

  • I think students are understanding now that the role of the ‘ticks’ is to help identify difficulties. I need to do a better job next year at explaining this from the offset. It might help to put it at the top of every lab sheet… I’ll think about this…

  • Imogen’s focus group went really well with 15 students participating. I still need to carefully reflect on the particular issues that were raised. I feel again that some need to be addressed through a better communication on my part (for example students wanted to have a list of alternative resources, there is such a list on the class site and I thought I’d mentioned it sufficiently but I will need to do that better).

  • There are a large number of tutors now and with that comes the tricky task of ensuring they understand their role. Jason was mentioned in Imogen’s focus group as doing a great job as “he didn’t just give the answer but pushed students to identify their difficulties”. I need to think very carefully about how to get this across to all the tutors (I think this has been addressed for the current year).

  • I’ve been thinking about the resources and am thinking that I might create a second set of videos. This would be a large quantity of work (3/4 weekends probably) but I think I could really help student further. In a way it seems logical to do after having run the course twice. This would further ‘flip the class’…

  • On another good note not entirely irrelevant to the class: code club is going well! The last session was a very busy one and students are actively participating which is great. Some first years are starting to work on the Euler problems and other formed a bit of a revision group for the upcoming class test they have… You can see the site (that has been really nicely stylised by the students):

October 25, 2014 12:00 AM

October 17, 2014

William Stein

A Non-technical Overview of the SageMathCloud Project

What problems is the SageMathCloud project trying to solve? What pain points does it address? Who are the competitors and what is the state of the technology right now?

What problems you’re trying to solve and why are these a problem?

  • Computational Education: How can I teach a course that involves mathematical computation and programming?
  • Computational Research: How can I carry out collaborative computational research projects?
  • Cloud computing: How can I get easy user-friendly collaborative access to a remote Linux server?

What are the pain points of the status quo and who feels the pain?

  • Student/Teacher pain:
    • Getting students to install software needed for a course on their computers is a major pain; sometimes it is just impossible, due to no major math software (not even Sage) supporting all recent versions of Windows/Linux/OS X/iOS/Android.
    • Getting university computer labs to install the software you need for a course is frustrating and expensive (time and money).
    • Even if computer labs worked, they are often being used by another course, stuffy, and students can't possibly do all their homework there, so computation gets short shrift. Lab keyboards, hardware, etc., all hard to get used to. Crappy monitors.
    • Painful confusing problems copying files around between teachers and students.
    • Helping a student or collaborator with their specific problem is very hard without physical access to their computer.
  • Researcher pain:
    • Making backups every few minutes of the complete state of everything when doing research often hard and distracting, but important for reproducibility.
    • Copying around documents, emailing or pushing/pulling them to revision control is frustrating and confusing.
    • Installing obscuring software is frustarting and distracting from the research they really want to do.
  • Everybody:
    • It is frustrating not having LIVE working access to your files wherever you are. (Dropbox/Github doesn't solve this, since files are static.)
    • It is difficult to leave computations running remotely.

Why is your technology poised to succeed?

  • When it works, SageMathCloud solves every pain point listed above.
  • The timing is right, due to massive improvements in web browsers during the last 3 years.
  • I am on full sabbatical this year, so at least success isn't totally impossible due to not working on the project.
  • I have been solving the above problems in less scalable ways for myself, colleagues and students since the 1990s.
  • SageMathCloud has many users that provide valuable feedback.
  • We have already solved difficult problems since I started this project in Summer 2012 (and launched first version in April 2013).

Who are your competitors?

There are no competitors with a similar range of functionality. However, there are many webapps that have some overlap in capabilities:
  • Mathematical overlap: Online Mathematica: "Bring Mathematica to life in the cloud"
  • Python overlap: Wakari: "Web-based Python Data Analysis"; also PythonAnywhere
  • Latex overlap: ShareLaTeX, WriteLaTeX
  • Web-based IDE's/terminals: target writing webapps (not research or math education):,,,
  • Homework: WebAssign and WebWork
Right now, SageMathCloud gives away for free far more than any other similar site, due to very substantial temporary financial support from Google, the NSF and others.

What’s the total addressable market?

Though our primary focus is the college mathematics classroom, there is a larger market:
  • Students: all undergrad/high school students in the world taking a course involving programming or mathematics
  • Teachers: all teachers of such courses
  • Researchers: anybody working in areas that involve programming or data analysis
Moreover, the web-based platform for computing that we're building lends itself to many other collaborative applications.

What stage is your technology at?

  • The site is up and running and has 28,413 monthly active users
  • There are still many bugs
  • I have a precise todo list that would take me at least 2 months fulltime to finish.

Is your solution technically feasible at this point?

  • Yes. It is only a matter of time until the software is very polished.
  • Morever, we have compute resources to support significantly more users.
  • But without money (from paying customers or investment), if growth continues at the current rate then we will have to clamp down on free quotas for users.

What technical milestones remain?

  • Infrastructure for creating automatically-graded homework problems.
  • Fill in lots of details and polish.

Do you have external credibility with technical/business experts and customers?

  • Business experts: I don't even know any business experts.
  • Technical experts: I founded the Sage math software, which is 10 years old and relatively well known by mathematicians.
  • Customers: We have no customers; we haven't offered anything for sale.

Market research?

  • I know about math software and its users as a result of founding the Sage open source math software project, NSF-funded projects I've been involved in, etc.

Is the intellectual property around your technology protected?

  • The IP is software.
  • The website software is mostly new Javascript code we wrote that is copyright Univ. of Washington and mostly not open source; it depends on various open source libraries and components.
  • The Sage math software is entirely open source.

Who are the team members to move this technology forward?

I am the only person working on this project fulltime right now.
  • Everything: William Stein -- UW professor
  • Browser client code: Jon Lee, Andy Huchala, Nicholas Ruhland -- UW undergrads
  • Web design, analytics: Harald Schilly -- Austrian grad student
  • Hardware: Keith Clawson

Why are you the ideal team?

  • We are not the ideal team.
  • If I had money maybe I could build the ideal team, leveraging my experience and connections from the Sage project...

by William Stein ( at October 17, 2014 01:04 PM

October 16, 2014

William Stein

Public Sharing in SageMathCloud, Finally

SageMathCloud (SMC) is a free (NSF, Google and UW supported) website that lets you collaboratively work with Sage worksheets, IPython notebooks, LaTeX documents and much, much more. All work is snapshotted every few minutes, and copied out to several data centers, so if something goes wrong with a project running on one machine (right before your lecture begins or homework assignment is due), it will pop up on another machine. We designed the backend architecture from the ground up to be very horizontally scalable and have no single points of failure.

This post is about an important new feature: You can now mark a folder or file so that all other users can view it, and very easily copy it to their own project.

This solves problems:
  • Problem: You create a "template" project, e.g., with pre-installed software, worksheets, IPython notebooks, etc., and want other users to easily be able to clone it as a new project. Solution: Mark the home directory of the project public, and share the link widely.

  • Problem: You create a syllabus for a course, an assignment, a worksheet full of 3d images, etc., that you want to share with a group of students. Solution: Make the syllabus or worksheet public, and share the link with your students via an email and on the course website. (Note: You can also use a course document to share files with all students privately.) For example...

  • Problem: You run into a problem using SMC and want help. Solution: Make the worksheet or code that isn't working public, and post a link in a forum asking for help.
  • Problem: You write a blog post explaining how to solve a problem and write related code in an SMC worksheet, which you want your readers to see. Solution: Make that code public and post a link in your blog post.
Here's a screencast.

Each SMC project has its own local "project server", which takes some time to start up, and serves files, coordinates Sage, terminal, and IPython sessions, etc. Public sharing completely avoids having anything to do with the project server -- it works fine even if the project server is not running -- it's always fast and there is no startup time if the project server isn't running. Moreover, public sharing reads the live files from your project, so you can update the files in a public shared directory, add new files, etc., and users will see these changes (when they refresh, since it's not automatic).
As an example, here is the cloud-examples github repo as a share. If you click on it (and have a SageMathCloud account), you'll see this:

What Next?

There is an enormous amount of natural additional functionality to build on top of public sharing.

For example, not all document types can be previewed in read-only mode right now; in particular, IPython notebooks, task lists, LaTeX documents, images, and PDF files must be copied from the public share to another project before people can view them. It is better to release a first usable version of public sharing before systematically going through and implementing the additional features needed to support all of the above. You can make complicated Sage worksheets with embedded images and 3d graphics, and those can be previewed before copying them to a project.
Right now, the only way to visit a public share is to paste the URL into a browser tab and load it. Soon the projects page will be re-organized so you can search for publicly shared paths, see all public shares that you have previously visited, who shared them, how many +1's they've received, comments, etc.

Also, I plan to eventually make it so public shares will be visible to people who have not logged in, and when viewing a publicly shared file or directory, there will be an option to start it running in a limited project, which will vanish from existence after a period of inactivity (say).

There are also dozens of details that are not yet implemented. For example, it would be nice to be able to directly download files (and directories!) to your computer from a public share. And it's also natural to share a folder or file with a specific list of people, rather than sharing it publicly. If somebody is viewing a public file and you change it, they should likely see the update automatically. Right now when viewing a share, you don't even know who shared it, and if you open a worksheet it can automatically execute Javascript, which is potentially unsafe.  Once public content is easily found, if somebody posts "evil" content publicly, there needs to be an easy way for users to report it.

Sharing will permeate everything

Sharing has been thought about a great deal during the last few years in the context of sites such as Github, Facebook, Google+ and Twitter. With SMC, we've developed a foundation for interactive collaborative computing in a browser, and will introduce sharing on top of that in a way that is motivated by your problems. For example, as with Github or Google+, when somebody makes a copy of your publicly shared folder, this copy should be listed (under "copies") and it could start out public by default. There is much to do.

One reason it took so long to release the first version of public sharing is that I kept imagining that sharing would happen at the level of complete projects, just like sharing in Github. However, when thinking through your problems, it makes way more sense in SMC to share individual directories and files. Technically, sharing at this level works works well for read-only access, not for read-write access, since projects are mapped to Linux accounts. Another reason I have been very hesitant to support sharing is that I've had enormous trouble over the years with spammers posting content that gets me in trouble (with my University -- it is illegal for UW to host advertisements). However, by not letting search engines index content, the motivation for spammers to post nasty content is greatly reduced.

Imagine publicly sharing recipes for automatically gradable homework problems, which use the full power of everything installed in SMC, get forked, improved, used, etc.

by William Stein ( at October 16, 2014 02:29 PM

October 13, 2014

Vince Knight

Reflecting on a first week of learning

This academic year Imogen Dunne (a final year student here at Cardiff University) is doing her final year project looking at evaluating student attitudes in my “Computing for Mathematics” course.

The idea is to use questionnaires, focus groups and interviews to evaluate (longitudinally) mathematics students’ attitudes towards:

  • Learning to code;
  • Learning some early mathematics;
  • Learning in a flipped class environment.

This has started off well, and last week Imogen kicked off a meeting with me by giving me homework. Something along the lines of:

I’d like to see the attitudes of everyone involved from every point of view. Could you perhaps write a diary/log at certain points during the year, describing how you feel things are going?

I was delighted with this idea and thought that I might as well blog these reflections, so here it goes:

A good first week.

The first thing I will say is I think this first week went really well. Students turned up to their lab sessions having almost universally carried out all the required work which was awesome. There was one major change with last year shifting some of the content from each week to the next. This had the effect of making the first week a bit lighter which I think has been a good thing.

Some labs could have been a bit ‘noisier’.

I want labs sessions to be noisy spaces with students talking to each other and figuring things out. This happened in most of the lab sessions I got to see but there was one or two where it felt more like a high school class from when I was still in high school: students looking at me as if I was a teacher evaluating whatever they were saying. I’m using 2nd year students as tutors this year (which I’m super excited about: more about that later) and perhaps I need to do a better job explaining exactly what it is that I want the labs to be. I think this has now been addressed.

The best flipped class meeting I’ve ever had.

I strive for a student centred learning environment. This isn’t always easy to obtain but the first class meeting went extremely well. I came in to the meeting expecting us to talk about the index method on strings (which finds the location of the first occurrence of a substring in a string) as this was the main piece of feedback I had from the labs as to the difficulties of the students. We started talking about it that but then the students really took over and wanted to know how to calculate the location of all the subtrings. This was a really awesome session as we went in to a particular thing in much more detail than we would have otherwise. Most importantly students would never have been able to have that conversation with me if they knew nothing about strings and the index function…

Not leaving anyone behind.

When teaching in a classic lecture based setting it’s extremely difficult to gain an understanding of how your students are doing. This flipped environment is all about finding out how students are doing, I am constantly getting feedback as to what students are having difficulty with. I have to make sure that students understand that that is what the labs are for, I’m constantly saying this but will continue to do so.

Further more there are some students who are having difficulties with the content, this is completely expected but the great thing about teaching using this approach is that I’ve already been able to identify them and will be meeting with them during my office hours to help (I’m always happy to help students who want to work).

Big thanks to the tutors.

Finally, this week of class would not have gone so smoothly if it weren’t for the great tutors who have helped in the labs. These include some of my colleagues who have gratefully given up their time, Jason Young who has done an awesome job organising the undergraduate tutors and most importantly the undergraduate tutors themselves. They all did a great job and I hope will also continue to learn and enjoy writing code.

This has been written down in a slight rush before my next set of labs but hopefully will be useful to Imogen (and indeed perhaps my students).

October 13, 2014 12:00 AM

October 12, 2014

Vince Knight

A playlist of 21 videos introducing LaTeX using SageMathCloud

About a year ago I put together 21 videos showing very basic LaTeX syntax. I used writelatex as the platform for those videos.

You can see that playlist of videos here.

This year I’m planning on using SageMathCloud when I teach Sage to my first years so I thought it was probably worth showing the students LaTeX in the same environment. So I’ve just finished redo’ing the same playlist of 21 videos but this time using SageMathCloud.

You can find that playlist here:

I’ve put all the tex files I created in the video in a github repository so they can be easily cloned in a SageMathCloud project using the https url:

October 12, 2014 12:00 AM

October 04, 2014

Vince Knight

A list of podcasts I listen to

I’ve recently gotten in to podcasts again (again) and really enjoy listening to most of these in the background of whatever it is I’m doing.

The podcasts I listen to probably fall in to the following categories:

  • Technology
  • Science
  • Sport

I’m sure I’m missing a bunch so I thought I’d write a blog post in the hope that kind readers would let me know what I should be listening to.


  • Tech News Today: One of many twit shows I like. I usually listen to this first thing in the morning.
  • This week in tech: I really like Leo Laporte and enjoy listening to this show as background noise to whatever it is that I’m doing.
  • Linux action show: I sometime find these guys ‘defend’ linux on the desktop a bit too fervently but I never fail to learn something new.
  • MacBreak Weekly: This has the same kind of vibe as ‘this week in tech’.
  • All about android: Again, this one feels like hanging out with friends (the probably with real friends being that I can’t work at the same time but I’m working on that).
  • Security Now: Another twit podcast, this one is about security and although I miss a lot about what is talked about I still enjoy it.
  • Le Rendez-Vous Tech: A French podcast very similar to this week in tech. Probably the only bit of exercise my French still gets.
  • Healthy Hacker: This is probably one of my favourite podcasts. Chris talks about code but also general fitness stuff.
  • Not really a podcast (low frequency) but it’s on my podcast catcher so I thought I’d list it here.


  • Freakonomics Radio: A great podcast, extremely well produced and often covering really interesting subjects (this one I usually try not to listen to in the background like most of the above).
  • The infinite monkey cage: Another really good job. Often very funny as well as informative.
  • This week in Science: this is always a nice listen to as they catch up on scientific stories that happened throughout the week.
  • Pythagoras Trousers: Cardiff University graduate Rhys Phillips does an excellent job with this podcast.
  • Math Ed Podcast: A really nice show that in every episode I’ve listened to was a really interesting interview of a mathematical education expert.
  • Stuff you should know: Not too sure if I should put this in this category but it’s a cool show where they explain a bunch of topics.


  • First Take: I don’t get to follow much US sports but I do enjoy listening to Slip and Stephen A. yell at each other.
  • Pardon the Interruption: Same as above really (although with less yelling).
  • Fighting Talk: This is a fun BBC show where sports pundits get points for being funny talking about stuff that happened in the past week.
  • Around the NFL: I like the NFL and when I have time to listen to this one I get to find out what’s happening.
  • Egg chasers podcast: A really great weekly roundup of rugby.
  • Scrum V Radio: A specific podcast about Welsh rugby.

One I’m forgetting in all of that is Relatively Prime which is a mathematics podcast that has just started a kickstarter campaign for a second season. I’ve never listened to it but though it could be worth mentioning.

Any good ones I’m missing?

October 04, 2014 12:00 AM

October 01, 2014

William Stein

SageMathCloud Course Management

by William Stein ( at October 01, 2014 01:05 PM

September 27, 2014

Sébastien Labbé

Abelian complexity of the Oldenburger sequence

The Oldenburger infinite sequence [O39] \[ K = 1221121221221121122121121221121121221221\ldots \] also known under the name of Kolakoski, is equal to its exponent trajectory. The exponent trajectory \(\Delta\) can be obtained by counting the lengths of blocks of consecutive and equal letters: \[ K = 1^12^21^22^11^12^21^12^21^22^11^22^21^12^11^22^11^12^21^22^11^22^11^12^21^12^21^22^11^12^21^12^11^22^11^22^21^12^21^2\ldots \] The sequence of exponents above gives the exponent trajectory of the Oldenburger sequence: \[ \Delta = 12211212212211211221211212\ldots \] which is equal to the original sequence \(K\). You can define this sequence in Sage:

sage: K = words.KolakoskiWord()
sage: K
word: 1221121221221121122121121221121121221221...
sage:          # delta returns the exponent trajectory
word: 1221121221221121122121121221121121221221...

There are a lot of open problem related to basic properties of that sequence. For example, we do not know if that sequence is recurrent, that is, all finite subword or factor (finite block of consecutive letters) always reappear. Also, it is still open to prove whether the density of 1 in that sequence is equal to \(1/2\).

In this blog post, I do some computations on its abelian complexity \(p_{ab}(n)\) defined as the number of distinct abelian vectors of subwords of length \(n\) in the sequence. The abelian vector \(\vec{w}\) of a word \(w\) counts the number of occurences of each letter: \[ w = 12211212212 \quad \mapsto \quad 1^5 2^7 \text{, abelianized} \quad \mapsto \quad \vec{w} = (5, 7) \text{, the abelian vector of } w \]

Here are the abelian vectors of subwords of length 10 and 20 in the prefix of length 100 of the Oldenburger sequence. The functions abelian_vectors and abelian_complexity are not in Sage as of now. Code is available at trac #17058 to be merged in Sage soon:

sage: prefix = words.KolakoskiWord()[:100]
sage: prefix.abelian_vectors(10)
{(4, 6), (5, 5), (6, 4)}
sage: prefix.abelian_vectors(20)
{(8, 12), (9, 11), (10, 10), (11, 9), (12, 8)}

Therefore, the prefix of length 100 has 3 vectors of subwords of length 10 and 5 vectors of subwords of length 20:

sage: p100.abelian_complexity(10)
sage: p100.abelian_complexity(20)

I import the OldenburgerSequence from my optional spkg because it is faster than the implementation in Sage:

sage: from slabbe import KolakoskiWord as OldenburgerSequence
sage: Olden = OldenburgerSequence()

I count the number of abelian vectors of subwords of length 100 in the prefix of length \(2^{20}\) of the Oldenburger sequence:

sage: prefix = Olden[:2^20]
sage: %time prefix.abelian_vectors(100)
CPU times: user 3.48 s, sys: 66.9 ms, total: 3.54 s
Wall time: 3.56 s
{(47, 53), (48, 52), (49, 51), (50, 50), (51, 49), (52, 48), (53, 47)}

Number of abelian vectors of subwords of length less than 100 in the prefix of length \(2^{20}\) of the Oldenburger sequence:

sage: %time L100 = map(prefix.abelian_complexity, range(100))
CPU times: user 3min 20s, sys: 1.08 s, total: 3min 21s
Wall time: 3min 23s
sage: from collections import Counter
sage: Counter(L100)
Counter({5: 26, 6: 26, 4: 17, 7: 15, 3: 8, 8: 4, 2: 3, 1: 1})

Let's draw that:

sage: labels = ('Length of factors', 'Number of abelian vectors')
sage: title = 'Abelian Complexity of the prefix of length $2^{20}$ of Oldenburger sequence'
sage: list_plot(L100, color='green', plotjoined=True, axes_labels=labels, title=title)

It seems to grow something like \(\log(n)\). Let's now consider subwords of length \(2^n\) for \(0\leq n\leq 12\) in the same prefix of length \(2^{20}\):

sage: %time L20 = [(2^n, prefix.abelian_complexity(2^n)) for n in range(20)]
CPU times: user 41 s, sys: 239 ms, total: 41.2 s
Wall time: 41.5 s
sage: L20
[(1, 2), (2, 3), (4, 3), (8, 3), (16, 3), (32, 5), (64, 5), (128, 9),
(256, 9), (512, 13), (1024, 17), (2048, 22), (4096, 27), (8192, 40),
(16384, 46), (32768, 67), (65536, 81), (131072, 85), (262144, 90), (524288, 104)]

I now look at subwords of length \(2^n\) for \(0\leq n\leq 23\) in the longer prefix of length \(2^{24}\):

sage: prefix = Olden[:2^24]
sage: %time L24 = [(2^n, prefix.abelian_complexity(2^n)) for n in range(24)]
CPU times: user 20min 47s, sys: 13.5 s, total: 21min
Wall time: 20min 13s
sage: L24
[(1, 2), (2, 3), (4, 3), (8, 3), (16, 3), (32, 5), (64, 5), (128, 9), (256,
9), (512, 13), (1024, 17), (2048, 23), (4096, 33), (8192, 46), (16384, 58),
(32768, 74), (65536, 98), (131072, 134), (262144, 165), (524288, 229),
(1048576, 302), (2097152, 371), (4194304, 304), (8388608, 329)]

The next graph gather all of the above computations:

sage: G = Graphics()
sage: legend = 'in the prefix of length 2^{}'
sage: G += list_plot(L24, plotjoined=True, thickness=4, color='blue', legend_label=legend.format(24))
sage: G += list_plot(L20, plotjoined=True, thickness=4, color='red', legend_label=legend.format(20))
sage: G += list_plot(L100, plotjoined=True, thickness=4, color='green', legend_label=legend.format(20))
sage: labels = ('Length of factors', 'Number of abelian vectors')
sage: title = 'Abelian complexity of Oldenburger sequence'
sage:'semilogx', 2), axes_labels=labels, title=title)

A linear growth in the above graphics with logarithmic \(x\) abcisse would mean a growth in \(\log(n)\). After those experimentations, my hypothesis is that the abelian complexity of the Oldenburger sequence grows like \(\log(n)^2\).


[O39]Oldenburger, Rufus (1939). "Exponent trajectories in symbolic dynamics". Transactions of the American Mathematical Society 46: 453–466. doi:10.2307/1989933

by Sébastien Labbé at September 27, 2014 10:00 PM

September 24, 2014

Vince Knight

Grey scale network graph colorings in Sage

This is a quick post following a request for some Sage help that a colleague asked for. It’s based on a quick fix and I’m wondering if someone might come up with a better way of doing this that I’ve missed or if it’s worth actually raising a ticket to incorporate something like it in Sage.

So my colleague is writing a book on Graph theory and recently started taking a look at Sage’s capacity to handle Graph theory stuff and colorings in particular. The issue was that said colleague ideally wanted grey scale pictures of the colorings (I’m guessing this is due to the publisher or something - I didn’t ask).

The following creates a bespoke graph and plots a coloring (ensures that adjacent vertices have different colors):

sage: P = graphs.PetersenGraph()
sage: c = P.coloring()
sage: c
[[1, 3, 5, 9], [0, 2, 6], [4, 7, 8]]

Now to get that in to grey scale we could of course open up inkscape or something similar and convert it but it would be nice to be able to directly use something like the matplotlib grey scale color map. This is in fact what I started to look for but with no success so I then started to look for how one converts an RGB tuple (3 floats corresponding to the makeup of a color) to something on a grey scale.

Turns out (see this stackoverflow question which leads to this wiki page), that for \(\text{rgb}=(r,g,b)\), the corresponding grey scale color is given by \(\text{grey}=(Y,Y,Y)\) where \(Y\) is given by:

where \(y\) is given by:

I genuinely have no understanding what so ever as to what that does but the idea is to make use of the Sage rainbow function which returns a given number of colors (very useful for creating plots when you don’t necessarily want to come up with all the color names).

sage: rainbow(10, 'rgbtuple')
[(1.0, 0.0, 0.0),
 (1.0, 0.6000000000000001, 0.0),
 (0.7999999999999998, 1.0, 0.0),
 (0.20000000000000018, 1.0, 0.0),
 (0.0, 1.0, 0.40000000000000036),
 (0.0, 1.0, 1.0),
 (0.0, 0.40000000000000036, 1.0),
 (0.1999999999999993, 0.0, 1.0),
 (0.8000000000000007, 0.0, 1.0),
 (1.0, 0.0, 0.5999999999999996)]

So here’s a function that takes the output of rainbow and maps it to grey scale:

def grey_rainbow(n, black=False):
    Return n greyscale colors
    if black:
        clrs = [0.299*clr[0] + 0.587 * clr[1] + 0.114 * clr[2] for clr in rainbow(n-2,'rgbtuple')]
        clrs = [0.299*clr[0] + 0.587 * clr[1] + 0.114 * clr[2] for clr in rainbow(n-1,'rgbtuple')]
    output = ['white']
    for c in clrs:
        if c <= 0.0031308:
            rgb = 12.92 * c
            rgb = (1.055*c^(1/2.4)-0.055)
    if black:
    return output

Note that we’re including the option to use black as one of the colours or not (it covers up the vertex labels on the corresponding plot if we do). We can then use that function to create our own partition coloring:

def grey_coloring(G, black=False):
    chromatic_nbr = G.chromatic_number()
    coloring = G.coloring()
    grey_colors = grey_rainbow(chromatic_nbr, black)
    d = {}
    for i, c in enumerate(grey_colors):
        d[c] = coloring[i]
    return P.graphplot(vertex_colors=d)

Here is how we can simply use the above to get a grey scale coloring of a graph:

P = graphs.PetersenGraph()
p = grey_coloring(P)

P = graphs.PetersenGraph()
p = grey_coloring(P,black=True)

Now what would be really nice would be to be able to just use any matplotlib color map in the Graph coloring. This might actually already be possible, I’ll fish through the Sage source code at some point and take a look (the awesome thing about Sage is that I can do that). Otherwise, it might just be a quick fix (and hopefully a less hacky one then above - I still laugh at the formulae I use and that seems to work), who nows I might even see if it’s worth opening an actual ticket for this.

Here is a Sage script with the above code

September 24, 2014 12:00 AM

September 21, 2014

Vince Knight

My thoughts on plotly (the github for plots)

A while ago I saw plotly appear on my G+ stream. People seemed excited about it but I was too busy to look at it properly and just assumed: must be some sort of new matplotlib: ain’t nobody got time for that!

Then, one of the guys from plotly reached out saying I should take a look. I took a brief glance and realised that this was nothing like a new matplotlib and in fact looked pretty cool. So I dutifully put it on my to do list but very much near the bottom.

I’m writing this sat in between sessions at PyconUK 2014. One of the talks on the first day was by Chris from plotly. He gave a great talk (which once the video link is up I’ll share here) and I immediately threw ‘check out plotly’ to the top of my to do list.

How I got started

import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import numpy as np
import plotly.plotly as py

n = 50
x, y, z, s, ew = np.random.rand(5, n)
c, ec = np.random.rand(2, n, 4)
area_scale, width_scale = 500, 5

fig, ax = plt.subplots()
sc = ax.scatter(x, y, c=c,

plot_url = py.plot_mpl(fig)

That code automatically creates the following plotly plot (which you can edit, zoom in etc…):

By ‘automatically’ I mean: ‘opens up web browser and your plots is there’!!

Doing something of my own

In my previous post I wrote about how to use Markov Chains to obtain the expected wait in a tandem qeue. Here’s a plot I put together that compared the analytical values to simulated values:

The code to obtain that particular plot is below:

# Libraries
from matplotlib.pyplot import plt
import csv
import plotly.plotly as py

# Get parameters

c1   = 12
N    = 9
c2   = 12
mu_1 = 1
mu_2 = .2
p    = .5

# Read analytical data
analytical_data = [[float(k) for k in row] for row in csv.reader(open('analytical.csv', 'r'))]

# Read simulation data
simulation_data = [[float(k) for k in row] for row in csv.reader(open('simulated.csv', 'r'))]

# Create the plot

fig = plt.figure()
ax = plt.subplot(111)
x_sim = [row[0] for row in simulation_data[::10]]  # The datasets have more data than I want to plot so skipping some values
y_sim = [row[1:] for row in simulation_data[::10]]
ax.boxplot(y_sim, positions=x_sim)
x_ana = [row[0] for row in analytical_data if row[0] <= max(x_sim)]
y_ana = [row[1] for row in analytical_data[:len(x_ana)]]
plt.xticks(range(0,int(max(x_ana) + max(int(max(x_ana)) / 10,1)), max(int(max(x_ana)) / 10,1)))
ax.set_ylabel('Mean expected wait')
title="$c_1=%s,\; N=%s,\; c_2=%s,\;\mu_1=%s,\; \mu_2=%s,\; p=%s $" % (c1, N, c1, mu_1, mu_2, p)

# Save the plot as a pdf file locally
plt.savefig("%s-%s-%s-%s-%s-%s.pdf" % (int(c1), int(N), int(c1), mu_1, mu_2, p), bbox_inches='tight')

plot_url = py.plot_mpl(fig)

If you’d like to repeat the above you can download the analytical and simulated datafiles.

The result of that can be seen here:

Further more that is just a ‘thing’ on my plotly profile so you can see it at this url:

Getting other formats

On that page I can tweak the graph if I wanted to and finally I can just grab the plot in whatever format I want by simply adding the correct format extension to the url:

My overall thoughts

So right now I’m just kind of excited about the possibilities (too many ideas to coherently filter out the good ones), there are also packages for R so I might try and get my students to play around with it in R when I teach it…

As a research tool, I think this will also be nice (it’s certainly the way to go). Although recently, I’ve been working remotely with two students and being able to throw a png of a plot in a hangout chat is pretty cool (and mobile friendly). So maybe that’s something the plotly guys could think about…

At the end of the day: this is an awesome tool. Plotly ‘abstractifies’ plots so that people using different packages/languages can still talk to each other. One of the big things I’m forgetting to talk about in detail is that there’s a web tool that allows you to change colors, change titles, mess with the data etc. That’s also a very cool collaborative tool obviously as I can imagine throwing up a plot that a co-author who doesn’t like code could then tweak.

Similarly (if/when) publications start using smarter formats (than things that are restricted by the need to be printed on paper) you could even just embed the plots like I’ve done here (so people could zoom, grab the data etc…). Here’s another way I could put that:

Papers are where plots go to die, they can go to plotly to live…

Woops, I’ve started blurting out some ideas… Hopefully they’re good ones.

I look forward to playing around with this tool some more (I need to see how it behaves with a Sage plot…).

September 21, 2014 12:00 AM

September 19, 2014

Vince Knight

Calculating the expected wait in a tandem queue with blocking using Sage

In this post I’ll describe a particular mathematical model that I’ve been working on for the purpose of a research paper. This is actually part of some work that I’ve done with +James Campbell an undergraduate who worked with me over the Summer.

Consider the system shown in the picture below:

This is a system composed of ‘two queues in tandem’ each with a number of servers \(c\) and a service rate \(\mu\). It is assumed here (as is often the case in queueing theory) that the service rate is exponentially distributed (in effect that the service time is random with mean \(1/\mu\)).

Individuals arrive at the queue with mean arrival rate \(\Lambda\) (once again exponentially distributed).

There is room for up to \(N\) individuals to wait for a free server at the first station. After their service at the first station is complete, individuals leave the system with probability \(p\), but if they don’t and there is no free place in the next station (ie there are not less than \(c_2\) individuals in the second service center) then they become blocked.

There are a vast array of applications of queueing systems like the above (in particular in the paper we’re working on we’re using it to model a healthcare system).

In this post I’ll describe how to use Markov chains to be able to describe the system and in particular how to get the expected wait for an arrival at the queue.

I have discussed Markov chains before (mainly posts on my old blog) and so I won’t go in to too much detail about the underlying theory (I think this is perhaps a slightly technical post so it is mainly aimed at people familiar with queueing theory but by all means ask in the comments if I can clarify anything).

The state space

One has to think about this carefully as it’s important to keep track not only of where individuals are but whether or not they are blocked. Based on that one might think of using a 3 dimensional state space: \((i,j,k)\) where \(i\) would denote the number of people at the first station, \(j\) the number at the second and \(k\) the number at the first who are blocked. This wouldn’t be a terrible idea but we can actually do things in a slightly neater and more compact way:

In the above \(i\) denotes the number of individuals in service or waiting for service at the first station and \(j\) denotes the number of individuals in service at the second station or blocked at the first station.

The continuous time transition rates between two states \((i_1,j_1)\) and \((i_2, j_2)\) are given by:

where \(\delta=(i_2,j_2)-(i_1,j_1)\).

Here’s a picture of the Markov Chain for \(N=c_1=c_2=2\):

The steady state probabilities

Using the above we can index our states and keep all the information in a matrix \(Q\), to obtain the steady state distributions of the chain (the probabilities of finding the queue in a given state) we then simply solve the following equation:

subject to \(\sum \pi = 1\).

Here’s how we can do this using Sage:

class Tandem_Queue():
        A class for an instance of the tandem_queue
        def __init__(self, c_1, N, c_2, Lambda, mu_1, mu_2, p):
            self.c_1 = c_1
            self.c_2 = c_2
            self.N = N
            self.Lambda = Lambda
            self.mu_1 = mu_1
            self.mu_2 = mu_2
            self.p = p
            self.m = c_1 + c_2 + 1
            self.n = c_2 + 1
            self.state_space = [(i, j)  for j in range(c_1 + c_2 + 1) for i in range(self.c_1 + self.N - max(j - self.c_2, 0) + 1)]
            if p == 1:  # Reduces state space in particular case of p = 1
                self.state_space = [state for state in self.state_space if state[1] == 0]
            Q = [[self.q(state1, state2) for state2 in self.state_space] for state1 in self.state_space]
            for i in range(len(Q)):
                Q[i][i] = - sum(Q[i])
            self.Q = matrix(QQ, Q)
            self.expected_wait_cache = {}

        def q(self, state1, state2):
            Returns the rate of transition between to given states.
            delta = list(vector(state2) - vector(state1))
            if delta == [1, 0]:
                return self.Lambda
            if delta == [-1, 1]:
                return min(self.c_1 - max(state1[1] - self.c_2, 0), state1[0]) * self.mu_1 * (1 - self.p)
            if delta == [-1, 0]:
                return min(self.c_1 - max(state1[1] - self.c_2, 0), state1[0]) * self.mu_1 * self.p
            if delta == [0, -1]:
                return min(state1[1], self.c_2) * self.mu_2
            return 0

        def pi(self):
            Solves linear system.
            A = transpose(self.Q).stack(vector([1 for state in self.state_space]))
            return A.solve_right(vector([0 for state in self.state_space] + [1]))

Most of the above is glorified book keeping but here’s a quick example showing what the above does and how it can be used. First let’s create an instance of our problem with \(N=c_1=c_2=2\), \(5\mu_2=2p=\mu_1=1\) and \(\Lambda=5\).

sage: small_example = tandem_queue(2,2,2,5,1,1/5,1/2)
sage: small_example.Q
22 x 22 dense matrix over Rational Field (use the '.str()' method to see the entries)

We see that if we return \(Q\) we get a 22 by 22 matrix which if you recall the picture above corresponds to the 22 states in that picture.

We can see the states by just typing:

sage: small_example.state_space
[(0, 0),
 (1, 0),
 (2, 0),
 (3, 0),
 (4, 0),
 (0, 1),
 (1, 1),
 (2, 1),
 (3, 1),
 (4, 1),
 (0, 2),
 (1, 2),
 (2, 2),
 (3, 2),
 (4, 2),
 (0, 3),
 (1, 3),
 (2, 3),
 (3, 3),
 (0, 4),
 (1, 4),
 (2, 4)]

If you check carefully they all correspond to the states of the picture above.

Now what we would like to know is the probability of being in any given state. To do this we need to solve the matrix equation \(\pi Q = 0\) such that \(\sum \pi=1\).

This is done in Sage (for any matrix Q) using the following:

sage: A = transpose(Q).stack(vector([1 for state in self.state_space]))
sage: A.solve_right(vector([0 for state in self.state_space] + [1]))

There we build a matrix A with an added column of 1s and then solve the corresponding equation using solve_right (note that we transposed the matrix). If you look at the class definition this was all defined earlier so we can in fact just run:

sage: small_example.pi()
(974420487508335328592/55207801002325145206717477, 6717141852060739142960/55207801002325145206717477, 25263720112088475982400/55207801002325145206717477, 107693117265184715581000/55207801002325145206717477, 499825193288571759140000/55207801002325145206717477, 7567657557556535357400/55207801002325145206717477, 50835142813671411162000/55207801002325145206717477, 177836071295654602905000/55207801002325145206717477, 638540135036394350075000/55207801002325145206717477, 2305924001256099701875000/55207801002325145206717477, 26439192416069771765000/55207801002325145206717477, 185599515623092483825000/55207801002325145206717477, 700026867396942548625000/55207801002325145206717477, 2256398553097737112500000/55207801002325145206717477, 4700830318953618984375000/55207801002325145206717477, 61385774570987050093750/55207801002325145206717477, 444444998037114715312500/55207801002325145206717477, 3393156381219452445312500/55207801002325145206717477, 15476151589322058007812500/55207801002325145206717477, 41152314633066177343750/55207801002325145206717477, 352285141439825390625000/55207801002325145206717477, 1826827211896183837890625/4246753923255780400516729)

That’s not super helpful displayed like that (all the arithmetic is being done exactly) so here’s a quick plot of the probabilities:

sage: p = list_plot(small_example.pi(), axes_labels = ['State', 'Probability'])
sage: p

That plot could be made a lot prettier an informative (by for example using the names of the states as xticks) but it will do for now. We can see from there for example that the most probable state of our queue (with the parameters we picked) is to be in the last state (see list above) which is \((2,4)\).

Out of interest here’s a plot when we change $\Lambda=1/2$ (a tenth of what we did above):

We see that now the most probable state is the sixth state (Python/Sage indexing starts at 0), which corresponds to \((0,1)\).

Here’s an animated plot of the steady state distribution for a larger system as \(\Lambda\) increases, displayed in a more informative way (the two dimensions corresponding to \(i\) and \(j\)):

All of that is very nice and interesting but where things get very useful is when trying to calculate the mean time one would expect to wait in a queue.

Mean expected wait in the queue

If we consider all states \((i,j)\in S\) only a subset of them will actually imply a wait:

  • If there are less than \(c_1\) individuals in the first station then anyone who arrives has direct access to a server;
  • If there are more than \(N+c_1\) individuals in the first station then anyone who arrives will be lost to the system.

With a little bit of thought (recalling what the \(i\)s and \(j\)s represent) we see that the states that incur a wait are given by:

Now if we know the expected wait when arriving in any state \((i,j) \in S_A\) we can get the mean wait as:

where \(w(i,j)\) denotes the expected time spent in any given state \((i,j)\). We are in essence, summing over arrivals that will have a wait and dividing by the probability of an individual not being lost to the system, which happens when \(i + \max(j - c_2, 0) < c_1 + N\).

Obtaining \(c(i,j)\) is relatively simple. We consider the ‘almost’ same Markov chain as above (except that we ignore arrivals). In this instance a jump from one step to another will only occur if a service occurs at the first station (with rate \(\min(c_1-\max(j-c_2,0),i)\mu_1\)) or if a service at the second station (with rate \(\min(c_2, j)\mu_2\)).

So the mean time spent in state \((i,j)\) is the inverse of the total exit rate:

Using that notion we are in effect discretizing the ‘ignore arrivals’ Markov chain. Once a transition occurs we can obtain the probability of where the service occurs:

  • The probability of the service being at the first station:
  • The probability of the service being at the second station:

We can use all of the above to put together a nice recursive formula for the expected wait \(w(i,j)\) in terms of the expected wait of states that are in effect getting closer and closer to having no wait:

where \(A\), is the set of states with no wait: \(i+\max(j-c_2,0) < c_1\)

Using the recursive formula is actually very easy to implement in code (we use a dictionary to cache calculated values to not make sure we don’t waste any time).

Here is a reduced version of the methods that need to be added to above to get this to work in Sage (you need to add self.expected_wait_cache = {} to the __init__ method):

def _pi_dict(self):
        Obtain a dictionary which indexes the states.
        self.pi_list = self.pi()
        return {state:self.pi_list[index] for index, state in enumerate(self.state_space)}

    def p_service_1(self, state):
        Returns the discretized probability of a service occurring at first station
        if self.p == 1:
            return 1
        return min(self.c_1 - max(state[1]- self.c_2, 0), state[0]) * self.mu_1 / (min(self.c_1 - max(state[1]- self.c_2, 0), state[0]) * self.mu_1 + min(self.c_2, state[1]) * self.mu_2)

    def p_service_2(self, state):
        Returns the discretized probability of a service occurring at second station
        if self.p == 1:
            return 0
        return  min(self.c_2, state[1]) * self.mu_2 / (min(self.c_1 - max(state[1]- self.c_2, 0), state[0]) * self.mu_1 + min(self.c_2, state[1]) * self.mu_2)

    def mean_time_in_state(self, state):
        Returns the mean time in any given state before a transition occurs
        return  1 / (min(self.c_1 - max(state[1]- self.c_2, 0), state[0]) * self.mu_1 + min(self.c_2, state[1]) * self.mu_2)

    def expected_wait(self, state):
        Function that returns the expected time till absorption for a given state
        if state in self.expected_wait_cache:
            return self.expected_wait_cache[state]
        if state not in self.state_space:  # If state outside of boundary. (Might not need this after new conditions below)
            return 0
        if state[0] + max(state[1] - self.c_2, 0) < self.c_1:  # If absorbing state
            self.expected_wait_cache[state] = 0
            return 0
        self.expected_wait_cache[state] =  (self.mean_time_in_state(state) + self.p * self.p_service_1(state) * self.expected_wait((state[0] - 1, state[1])))
        if (state[0] - 1, state[1] + 1) in self.state_space:
            self.expected_wait_cache[state] += (1-self.p) * self.p_service_1(state) * self.expected_wait((state[0] - 1, state[1] + 1))
        if (state[0], state[1] - 1) in self.state_space:
            self.expected_wait_cache[state] += self.p_service_2(state) * self.expected_wait((state[0], state[1] - 1))
        return self.expected_wait_cache[state]

    def mean_expected_wait(self):
        Returns the mean wait
        self.pi_dict = self._pi_dict()
        accepting_states = [state for state in [s for s in self.state_space if s[0] + max(s[1] - self.c_2, 0) < self.c_1 + self.N]]
        prob_of_accepting = sum([self.pi_dict[state] for state in accepting_states])
        return sum([self.expected_wait(state) * self.pi_dict[state] for state in accepting_states]) / prob_of_accepting

Using all of the above we can get the expected wait for our system:

sage: small_example = Tandem_Queue(2,2,2,1/2,1,1/5,1/2)
sage: small_example.mean_expected_wait()
sage: n(_)

Below is a plot showing the effect on the mean wait as demand increases in the plot below for a large system:

What that plot shows is the calculated values (solid blue) line going through box plots of the simulated value. Perhaps in another blog post some day I’ll write about how to simulate the above but I think that’s probability sufficient for now.

If it’s of interest all of the above code can be downloaded here or at this gist.

September 19, 2014 12:00 AM

August 27, 2014

Sébastien Labbé

slabbe-0.1.spkg released

These is a summary of the functionalities present in slabbe-0.1 optional Sage package. It depends on version 6.3 of Sage because it uses RecursivelyEnumeratedSet code that was merged in 6.3. It contains modules on digital geometry, combinatorics on words and more.

Install the optional spkg (depends on sage-6.3):

sage -i

In each of the example below, you first have to import the module once and for all:

sage: from slabbe import *

To construct the image below, make sure to use tikz package so that view is able to compile tikz code when called:

sage: latex.add_to_preamble("\\usepackage{tikz}")
sage: latex.extra_preamble()

Draw the part of a discrete plane

sage: p = DiscretePlane([1,pi,7], 1+pi+7, mu=0)
sage: d = DiscreteTube([-5,5],[-5,5])
sage: I = p & d
sage: I
Intersection of the following objects:
Set of points x in ZZ^3 satisfying: 0 <= (1, pi, 7) . x + 0 < pi + 8
DiscreteTube: Preimage of [-5, 5] x [-5, 5] by a 2 by 3 matrix
sage: clip = d.clip()
sage: tikz = I.tikz(clip=clip)
sage: view(tikz, tightpage=True)

Draw the part of a discrete line

sage: L = DiscreteLine([-2,3], 5)
sage: b = DiscreteBox([0,10], [0,10])
sage: I = L & b
sage: I
Intersection of the following objects:
Set of points x in ZZ^2 satisfying: 0 <= (-2, 3) . x + 0 < 5
[0, 10] x [0, 10]
sage: I.plot()

Double square tiles

This module was developped for the article on the combinatorial properties of double square tiles written with Ariane Garon and Alexandre Blondin Massé [BGL2012]. The original version of the code was written with Alexandre.

sage: D = DoubleSquare((34,21,34,21))
sage: D
Double Square Tile
  w0 = 3032321232303010303230301012101030   w4 = 1210103010121232121012123230323212
  w1 = 323030103032321232303                w5 = 101212321210103010121
  w2 = 2321210121232303232123230301030323   w6 = 0103032303010121010301012123212101
  w3 = 212323032321210121232                w7 = 030101210103032303010
(|w0|, |w1|, |w2|, |w3|) = (34, 21, 34, 21)
(d0, d1, d2, d3)         = (42, 68, 42, 68)
(n0, n1, n2, n3)         = (0, 0, 0, 0)
sage: D.plot()
sage: D.extend(0).extend(1).plot()

We have shown that using two invertible operations (called SWAP and TRIM), every double square tiles can be reduced to the unit square:

sage: D.plot_reduction()

The reduction operations are:

sage: D.reduction()
['SWAP_1', 'TRIM_1', 'TRIM_3', 'SWAP_1', 'TRIM_1', 'TRIM_3', 'TRIM_0', 'TRIM_2']

The result of the reduction is the unit square:

sage: unit_square = D.apply(D.reduction())
sage: unit_square
Double Square Tile
  w0 =     w4 =
  w1 = 3   w5 = 1
  w2 =     w6 =
  w3 = 2   w7 = 0
(|w0|, |w1|, |w2|, |w3|) = (0, 1, 0, 1)
(d0, d1, d2, d3)         = (2, 0, 2, 0)
(n0, n1, n2, n3)         = (0, NaN, 0, NaN)
sage: unit_square.plot()

Since SWAP and TRIM are invertible operations, we can recover every double square from the unit square:

sage: E = unit_square.extend(2).extend(0).extend(3).extend(1).swap(1).extend(3).extend(1).swap(1)
sage: D == E

Christoffel graphs

This module was developped for the article on a d-dimensional extension of Christoffel Words written with Christophe Reutenauer [LR2014].

sage: G = ChristoffelGraph((6,10,15))
sage: G
Christoffel set of edges for normal vector v=(6, 10, 15)
sage: tikz = G.tikz_kernel()
sage: view(tikz, tightpage=True)

Bispecial extension types

This module was developped for the article on the factor complexity of infinite sequences genereated by substitutions written with Valérie Berthé [BL2014].

The extension type of an ordinary bispecial factor:

sage: L = [(1,3), (2,3), (3,1), (3,2), (3,3)]
sage: E = ExtensionType1to1(L, alphabet=(1,2,3))
sage: E
  E(w)   1   2   3
   1             X
   2             X
   3     X   X   X
 m(w)=0, ordinary
sage: E.is_ordinaire()

Creation of a strong-weak pair of bispecial words from a neutral not ordinaire word:

sage: p23 = WordMorphism({1:[1,2,3],2:[2,3],3:[3]})
sage: e = ExtensionType1to1([(1,2),(2,3),(3,1),(3,2),(3,3)], [1,2,3])
sage: e
  E(w)   1   2   3
   1         X
   2             X
   3     X   X   X
 m(w)=0, not ord.
sage: A,B = e.apply(p23)
sage: A
  E(3w)   1   2   3
    2         X   X
    3     X   X   X
 m(w)=1, not ord.
sage: B
  E(23w)   1   2   3
    1          X
    3              X
 m(w)=-1, not ord.

Fast Kolakoski word

This module was written for fun. It uses cython implementation inspired from the 10 lines of C code written by Dominique Bernardi and shared during Sage Days 28 in Orsay, France, in January 2011.

sage: K = KolakoskiWord()
sage: K
word: 1221121221221121122121121221121121221221...
sage: %time K[10^5]
CPU times: user 1.56 ms, sys: 7 µs, total: 1.57 ms
Wall time: 1.57 ms
sage: %time K[10^6]
CPU times: user 15.8 ms, sys: 30 µs, total: 15.8 ms
Wall time: 15.9 ms
sage: %time K[10^8]
CPU times: user 1.58 s, sys: 2.28 ms, total: 1.58 s
Wall time: 1.59 s
sage: %time K[10^9]
CPU times: user 15.8 s, sys: 12.4 ms, total: 15.9 s
Wall time: 15.9 s

This is much faster than the Python implementation available in Sage:

sage: K = words.KolakoskiWord()
sage: %time K[10^5]
CPU times: user 779 ms, sys: 25.9 ms, total: 805 ms
Wall time: 802 ms


[BGL2012]A. Blondin Massé, A. Garon, S. Labbé, Combinatorial properties of double square tiles, Theoretical Computer Science 502 (2013) 98-117. doi:10.1016/j.tcs.2012.10.040
[LR2014]Labbé, Sébastien, and Christophe Reutenauer. A d-dimensional Extension of Christoffel Words. arXiv:1404.4021 (April 15, 2014).
[BL2014]V. Berthé, S. Labbé, Factor Complexity of S-adic sequences generated by the Arnoux-Rauzy-Poincaré Algorithm. arXiv:1404.4189 (April, 2014).

by Sébastien Labbé at August 27, 2014 04:53 PM

Releasing slabbe, my own Sage package

Since two years I wrote thousands of line of private code for my own research. Each module having between 500 and 2000 lines of code. The code which is the more clean corresponds to code written in conjunction with research articles. People who know me know that I systematically put docstrings and doctests in my code to facilitate reuse of the code by myself, but also in the idea of sharing it and eventually making it public.

I did not made that code into Sage because it was not mature enough. Also, when I tried to make a complete module go into Sage (see #13069 and #13346), then the monstrous never evolving #12224 became a dependency of the first and the second was unofficially reviewed asking me to split it into smaller chunks to make the review process easier. I never did it because I spent already too much time on it (making a module 100% doctested takes time). Also, the module was corresponding to a published article and I wanted to leave it just like that.

Getting new modules into Sage is hard

In general, the introduction of a complete new module into Sage is hard especially for beginners. Here are two examples I feel responsible for: #10519 is 4 years old and counting, the author has a new work and responsabilities; in #12996, the author was decouraged by the amount of work given by the reviewers. There is a lot of things a beginner has to consider to obtain a positive review. And even for a more advanced developper, other difficulties arise. Indeed, a module introduces a lot of new functions and it may also introduce a lot of new bugs... and Sage developpers are sometimes reluctant to give it a positive review. And if it finally gets a positive review, it is not available easily to normal users of Sage until the next release of Sage.

Releasing my own Sage package

Still I felt the need around me to make my code public. But how? There are people (a few of course but I know there are) who are interested in reproducing computations and images done in my articles. This is why I came to the idea of releasing my own Sage package containing my public research code. This way both developpers and colleagues that are user of Sage but not developpers will be able to install and use my code. This will make people more aware if there is something useful in a module for them. And if one day, somebody tells me: "this should be in Sage", then I will be able to say : "I agree! Do you want to review it?".

Old style Sage package vs New sytle git Sage package

Then I had to chose between the old and the new style for Sage packages. I did not like the new style, because

  • I wanted the history of my package to be independant of the history of Sage,
  • I wanted it to be as easy to install as sage -i slabbe,
  • I wanted it to work on any recent enough version of Sage,
  • I wanted to be able to release a new version, give it to a colleague who could install it right away without changing its own Sage (i.e., updating the checksums).

Therefore, I choose the old style. I based my work on other optional Sage packages, namely the SageManifolds spkg and the ore_algebra spkg.

Content of the initial version

The initial version of the slabbe Sage package has modules concerning four topics: Digital geometry, Combinatorics on words, Combinatorics and Python class inheritance.


For installation or for release notes of the initial version of the spkg, consult the slabbe spkg section of the Sage page of this website.

by Sébastien Labbé at August 27, 2014 04:48 PM

William Stein

What is SageMathCloud: let's clear some things up

[PDF version of this blog post]
"You will have to close source and commercialize Sage. It's inevitable." -- Michael Monagan, cofounder of Maple, told me this in 2006.
SageMathCloud (SMC) is a website that I first launched in April 2013, through which you can use Sage and all other open source math software online, edit Latex documents, IPython notebooks, Sage worksheets, track your todo items, and many other types of documents. You can write, compile, and run code in most programming languages, and use a color command line terminal. There is realtime collaboration on everything through shared projects, terminals, etc. Each project comes with a default quota of 5GB disk space and 8GB of RAM.

SMC is fun to use, pretty to look at, frequently backs up your work in many ways, is fault tolerant, encourages collaboration, and provides a web-based way to use standard command-line tools.

The Relationship with the SageMath Software

The goal of the SageMath software project, which I founded in 2005, is to create a viable free open source alternative to Magma, Mathematica, Maple, and Matlab. SMC is not mathematics software -- instead, SMC is best viewed by analogy as a browser-based version of a Linux desktop environment like KDE or Gnome. The vast majority of the code we write for SMC involves text editor issues (problems similar to those confronted by Emacs or Vim), personal information management, support for editing LaTeX documents, terminals, file management, etc. There is almost no mathematics involved at all.

That said, the main software I use is Sage, so of course support for Sage is a primary focus. SMC is a software environment that is being optimized for its users, who are mostly college students and teachers who use Sage (or Python) in conjunction with their courses. A big motivation for the existence of SMC is to make Sage much more accessible, since growth of Sage has stagnated since 2011, with the number one show-stopper obstruction being the difficulty of students installing Sage.

Sage is Failing

Measured by the mission statement, Sage has overall failed. The core goal is to provide similar functionality to Magma (and the other Ma's) across the board, and the Sage development model and community has failed to do this across the board, since after 9 years, based on our current progress, we will never get there. There are numerous core areas of research mathematics that I'm personally familiar with (in arithmetic geometry), where Sage has barely moved in years and Sage does only a few percent of what Magma does. Unless there is a viable plan for the areas to all be systematically addressed in a reasonable timeframe, not just with arithmetic geometry in Magma, but with everything in Mathematica, Maple., etc, we are definitely failing at the main goal I have for the Sage math software project.

I have absolutely no doubt that money combined with good planning and management would make it possible to achieve our mission statement. I've seen this hundreds of times over at a small scale at Sage Days workshops during the last decade. And let's not forget that with very substantial funding, Linux now provides a viable free open source alternative to Microsoft Windows. Just providing Sage developers with travel expenses (and 0 salary) is enough to get a huge amount done, when possible. But all my attempts with foundations and other clients to get any significant funding, at even the level of 1% of the funding that Mathematica gets each year, has failed. For the life of the Sage project, we've never got more than maybe 0.1% of what Mathematica gets in revenue. It's just a fact that the mathematics community provides Mathematica $50+ million a year, enough to fund over 600 fulltime positions, and they won't provide enough to fund one single Sage developer fulltime.

But the Sage mission statement remains, and even if everybody else in the world gives up on it, I HAVE NOT. SMC is my last ditch strategy to provide resources and visibility so we can succeed at this goal and give the world a viable free open source alternative to the Ma's. I wish I were writing interesting mathematical software, but I'm not, because I'm sucking it up and playing the long game.

The Users of SMC

During the last academic year (e.g., April 2014) there were about 20K "monthly active users" (as defined by Google Analytics), 6K weekly active users, and usually around 300 simultaneous connected users. The summer months have been slower, due to less teaching.

Numerically most users are undergraduate students in courses, who are asked to use SMC in conjunction with a course. There's also quite a bit of usage of SMC by people doing research in mathematics, statistics, economics, etc. -- pretty much all computational sciences. Very roughly, people create Sage worksheets, IPython notebooks, and Latex documents in somewhat equal proportions.

What SMC runs on

Technically, SMC is a multi-datacenter web application without specific dependencies on particular cloud provider functionality. In particular, we use the Cassandra database, and custom backend services written in Node.js (about 15,000 lines of backend code). We also use Amazon's Route 53 service for geographically aware DNS. There are two racks containing dedicated computers on opposites sides of campus at University of Washington with 19 total machines, each with about 1TB SSD, 4TB+ HDD, and 96GB RAM. We also have dozens of VM's running at 2 Google data centers to the east.

A substantial fraction of the work in implementing SMC has been in designing and implementing (and reimplementing many times, in response to real usage) a robust replicated backend infrastructure for projects, with regular snapshots and automatic failover across data centers. As I write this, users have created 66677 projects; each project is a self-contained Linux account whose files are replicated across several data centers.

The Source Code of SMC

The underlying source of SMC, both the backend server and frontend client, is mostly written in CoffeeScript. The frontend (which is nearly 20,000 lines of code) is implemented using the "progressive refinement" approach to HTML5/CSS/Javascript web development. We do not use any Javascript single page app frameworks, though we make heavy use of Bootstrap3 and jQuery. All of the library dependencies of SMC, e.g., CodeMirror, Bootstrap, jQuery, etc. for SMC are licensed under very permissive BSD/MIT, etc. libraries. In particular, absolutely nothing in the Javascript software stack is GPL or AGPL licensed. The plan is that any SMC source code that will be open sourced will be released under the BSD license. Some of the SMC source code is not publicly available, and is owned by University of Washington. But other code, e.g., the realtime sync code, is already available.
Some of the functionality of SMC, for example Sage worksheets, communicate with a separate process via a TCP connection. That separate process is in some cases a GPL'd program such as Sage, R, or Octave, so the viral nature of the GPL does not apply to SMC. Also, of course the virtual machines are running the Linux operating system, which is mostly GPL licensed. (There is absolutely no AGPL-licensed code anywhere in the picture.)

Note that since none of the SMC server and client code links (even at an interpreter level) with any GPL'd software, that code can be legally distributed under any license (e.g., from BSD to commercial).
Also we plan to create a fully open source version of the Sage worksheet server part of SMC for inclusion with Sage. This is not our top priority, since there are several absolutely critical tasks that still must be finished first on SMC, e.g., basic course management.

The SMC Business Model

The University of Washington Center for Commercialization (C4C) has been very involved and supportive since the start of the projects. There are no financial investors or separate company; instead, funding comes from UW, some unspent grant funds that were about to expire, and a substantial Google "Academic Education Grant" ($60K). Our first customer is the "US Army Engineer Research and Development Center", which just started a support/license agreement to run their own SMC internally. We don't currently offer a SaaS product for sale yet -- the options for what can be sold by UW are constrained, since UW is a not-for-profit state university. Currently users receive enhancements to their projects (e.g., increased RAM or disk space) in exchange for explaining to me the interesting research or teaching they are doing with SMC.

The longterm plan is to start a separate for-profit company if we build a sufficient customer base. If this company is successful, it would also support fulltime development of Sage (e.g., via teaching buyouts for faculty, support of students, etc.), similar to how Magma (and Mathematica, etc.) development is funded.

In conclusion, in response to Michael Monagan, you are wrong. And you are right.

by William Stein ( at August 27, 2014 07:55 AM

You don't really think that Sage has failed, do you?

I just received an email from a postdoc in Europe, and very longtime contributor to the Sage project.  He's asking for a letter of recommendation, since he has to leave the world of mathematical software development (after a decade of training and experience), so that he can take a job at hedge fund.  He ends his request with the question:

> P.S. You don't _really_ think that Sage has failed, do you?

After almost exactly 10 years of working on the Sage project, I absolutely do think it has failed to accomplish the stated goal of the mission statement: "Create a viable free open source alternative to Magma, Maple, Mathematica and Matlab.".     When it was only a few years into the project, it was really hard to evaluate progress against such a lofty mission statement.  However, after 10 years, it's clear to me that not only have we not got there, we are not going to ever get there before I retire.   And that's definitely a failure.   

Here's a very rough quote I overheard at lunch today at Sage Days 61, from John Voight, who wrote much quaternion algebra code in Magma: "I'm making a list of what is missing from Sage that Magma has for working with quaternion algebras.  However, it's so incredibly daunting, that I don't want to put in everything.  I've been working on Magma's quaternion algebras for over 10 years, as have several other people.  It's truly daunting how much functionality Magma has compared to Sage."

The only possible way Sage will not fail at the stated mission is if I can get several million dollars a year in money to support developers to work fulltime on implementing interesting core mathematical algorithms.  This is something that Magma has had for over 20 years, and that Maple, Matlab, and Mathematica also have.   That I don't have such funding is probably why you are about to take a job at a hedge fund.    If I had the money, I would try to hire a few of the absolute best people (rather than a bunch of amateurs), people like you, Robert Bradshaw, etc. -- we know who is good. (And clearly I mean serious salaries, not grad student wages!)

So yes, I think the current approach to Sage has failed.    I am going to try another approach, namely SageMathCloud.  If it works, maybe the world will get a free open source alternative to Magma, Mathematica, etc.  Otherwise, maybe the world never ever will.      If you care like I do about having such a thing, and you're teaching course, or whatever, maybe try using SageMathCloud.   If enough people use SageMathCloud for college teaching, then maybe a business model will emerge, and Sage will get proper funding.   

by William Stein ( at August 27, 2014 07:52 AM

Vince Knight

A Sneak Preview of Game Theory in Sage (2/3): Matching Games

In my previous post here I described some of the Sage development that +James Campbell and I spent a lot of time this Summer working on. In that post I described some work that has subsequently been accepted and included in the latest release of Sage (here’s the latest changlog): code to calculate the Shapley value.

In this post I’ll talk about the second of 3 tickets that James and I worked on: looking at Matching games. This has not actually been reviewed yet so please do help us get this code in to Sage by taking a look at the ticket: 16331.

What is a matching game?

One of the best explanations of a matching game (also called the stable marriage problem) can be found in this video. That video really is awesome but it might be a bit long (it’s 25 minutes) so this very short video I threw together for a class I teach might be of interest (it is no where near as good as the previous one but it’s 3 minutes long).

Basically a matching game attempts to create links between two groups of people (referred to as suitors and reviewers) in such a way as no one wants to break their link:

In the above picture we see the preferences of the suitors and the reviewers. So \(c\), prefers \(B\) to \(A\), and \(A\) to \(C\).

Here is the actual definition of a stable matching that I give my students:

A matching game of size \(N\) is defined by two disjoint sets \(S\) and \(R\) or suitors and reviewers of size \(N\). Associated to each element of \(S\) and \(R\) is a preference list:

A matching \(M\) is a any bijection between \(S\) and \(R\). If \(s\in S\) and \(r\in R\) are matched by \(M\) we denote:

The above image defines a matching game, one possible matching could be given below:

It’s immediate to note however that \(B\) and \(c\) prefer each other to their current matching: so the above matching is unstable. In that example \((B,c)\) is called a ‘blocking pair’.

Luckily Gale and Shapley obtained an algorithm that guarantees a stable matching and this is what James and I put together in to Sage.

First, let’s define a matching game:

sage: suitr_pref = {'a': ('B', 'A', 'C'),
....:               'b': ('B', 'C', 'A'),
....:               'c': ('A', 'B', 'C')}
sage: reviewr_pref = {'A': ('a', 'b', 'c'),
....:                 'B': ('a', 'c', 'b'),
....:                 'C': ('b', 'c', 'a')}
sage: m = MatchingGame([suitr_pref, reviewr_pref])

You can see that python dictionaries are used for the functions \(f\) and \(g\) described above (the suitor preferences).

If you tab complete after typing m. you can see some of the methods and attributes associated with the MatchingGame class:

sage: m.
m.add_reviewer  m.bi_partite    m.db            m.dumps         m.rename        m.reviewers     m.solve         m.version
m.add_suitor    m.category      m.dump          m.plot          m.reset_name          m.suitors

I won’t go in to much of the details of that year but you can get some help on anyone of those by typing ? after one of them (below you can see some of the output):

sage: m.solve?

Type:            instancemethod
File:            /Users/vince/sage/local/lib/python2.7/site-packages/sage/game_theory/
Definition:      m.solve(self, invert=False)
   Computes a stable matching for the game using the Gale-Shapley


Let’s give that method a spin (as you can see it’ll use the Gale-Shapley algorithm).

sage: m.solve()
{'C': ['b'], 'a': ['B'], 'b': ['C'], 'A': ['c'], 'c': ['A'], 'B': ['a']}

We see that a matching has been obtained. You can see the corresponding matching here:

Another nice method that we implemented is to use the awesome graph theory stuff that’s in Sage so you can obtain the corresponding bi-partite graph:

sage: p = m.bi_partite()
Bipartite graph on 6 vertices

You can see the corresponding plot here:

All of the above has not been reviewed yet so if you do have any comments they’d be very gratefully received. If you actually went over to trac and took a look at it there that would be great but otherwise just commenting here would be awesome.

This is the second in a series of 3 posts that I’ll get around to writing, in the next one I’ll cover ticket 16333: Normal Form Game. This is the biggest contribution by James as it involved interfacing with two other packages and also coding up a bespoke support enumeration algorithm.

August 27, 2014 12:00 AM