Meeting Julia

In my last post I mentioned the new language Julia. It deserves more than a single paragraph, so I thought I’d walk through a problem, and tackle it with the language.

The problem is a stock standard one: investigating the fitting of an SIR model of disease spread to an outbreak of influenza at an English boarding school. Of the 763 pupils, on the first day only one was infected, and the number of infected pupils changed over 14 days as:

1   2   3   4   5   6   7   8   9  10  11  12  13  14
3   6  25  73 222 294 258 237 191 125  69  27  11   4

We will try to find parameters \beta and \gamma so that when the equations

\displaystyle{\frac{dS}{dt}=-\beta IS}

\displaystyle{\frac{dI}{dt}=\beta IS-\gamma I}

\displaystyle{\frac{dR}{dt}=\gamma I}

are solved, the values of I fit the data.

Note that the sum of all the right hand sides is zero; this indicates that S+I+R is constant, which for this model is the population. (Other models allow for births, deaths and other population changes).

Finding the parameters can be done in several steps:

  1. Set up a function for the system of differential equations.
  2. From that function, create a second function which is the sum of squared differences between the data, and the computed values from the equations.
  3. Minimize that sum.

The first function can be created in Julia (following the example here) as:

function SIR(t,x,p)
     S = x[1]
     I = x[2]
     R = x[3]
     beta = p[1]
     gamma = p[2]
     [-beta*S*I, beta*S*I-gamma*I, gamma*I]
end

The default behaviour of functions in Julia is to return the line before the “end” statement; thus in a function such as this one, there is no need for a “return” statement (although one exists).

To solve a system of equations we need access to the ODE package, which I’ll assumed you’ve installed with

Pkg.add()

This can be loaded into a running Julia session either with

using ODE

which adds all the ODE functions into the top namespace, or as

import ODE

which adds all the functions in the ODE namespace. So to call the Dormand-Prince method, you would use

ode45

if all methods were in the top namespace, or as

ODE.ode45

in the second case.

So, for example, we could attempt to solve the equations as:

p = [0.01, 0.1]
t, x = ODE.ode45((t,x)->SIR(t,x,p),[0:0.1:14],[762.,1.,0.]);

To plot the functions, we’ll use the plotting library Winston (Winston, Julia from 1984 – geddit?), which is fairly basic, but enough for our simple needs:

import Winston
wn = Winston

wn.plot(t,x[:,1],"b",t,x[:,2],"g",t,x[:,3],"r",[1:14],data',"k*")

SIR_julia

The green curve corresponds to I, the number of infected individuals at a given time, and the asterisks to the data. Clearly these parameters do not provide a very good fit.

The next step is to create a sum of squares function to minimize. To do this, we will solve the ODE system with the ode4 method which uses a fixed step size. This means we can ensure that there are computed values for I at all the integer points:

function ss(b::Vector)
    data = [3 6 25 73 222 294 258 237 191 125 69 27 11 4];
    t,x = ode.ode4((t,x)->SIR(t,x,b),[0:0.05:14],[762.,1.,0.]);             
    sum((data-x[21:20:281,2]').^2)
    end

Now this function can be optimized, using a method from the Optim package:

import Optim

Optim.optimize(ss,[0.001,1],method=:cg)

The “method” variable here defaults to the Nelder-Mead method, but for this particular function I found that the conjugate-gradient method gave the best results.

This produces quite a bit of information, of which the salient bits are:

 * Minimum: .0021806376138654117
.44517197444683443

 * Value of Function at Minimum: 4507.078964

So, let’s try these parameters and see what happens:

p = [.0021806376138654117, .44517197444683443];
t, x = ODE.ode45((t,x)->SIR(t,x,p),[0:0.1:14],[762.,1.,0.]);
wn.plot(t,x[:,1],"b",t,x[:,2],"g",t,x[:,3],"r",[1:14],data',"k*")

And here’s the result:
SIR_julia2

As you can see, the fit is remarkably good.

As far as problems go, this is not particularly hard – conceptually speaking it’s probably at an undergraduate level. On the other hand, it’s not completely trivial, and does seem to me to give the computational environment a bit of a work-out. And Julia solves it with ease.

The best Matlab alternative (3)

Over two years ago I wrote The best Matlab alternative with a follow-up a bit later, which seem to have engendered a far amount of discussion.

Well, things have moved on in the computational world, and a user is now spoiled for choice.

Octave has just had version 3.8 released; this version comes with a GUI (which is not enabled by default). However, the GUI is a very fine thing indeed:

octave_gui

As you see, it has all the bells and whistles you may want: command window and history, variable browser, file browser, editor, and a documentation reader. The GUI will become standard as of version 4.0.

I think that Octave has come along in leaps and bounds over the past few years, and as a drop-in Matlab replacement (if that’s what you want) I don’t think it can be bettered.

Scilab is at version 5.4.1, and is a mature product. I have come across a few comments in newgroups to the effect that Scilab is better than Octave for some advanced numerical work. However, I’m not sure if these are just opinions or are backed up with solid data and timings. I used to like Scilab a great deal, but now I don’t see that it offers anything over Octave other than Xcos (a dynamic systems modeller, similar to Simulink), and a few niche extras, such as solutions of boundary value problems.

Python was a language I hadn’t considered in my first post, but with the addition of NumPy and SciPy (plus numerous add-on packages) becomes a definite contender for serious number-crunching. The web is full of stories about people who have dumped Matlab for Python. And of course with Python you get a lovely programming language, and extensions such as Cython, which gives you C-type speed. To get an idea of the packages, see here. Python has a interactive shell IPython which adds enormously to Python’s ease of use.

Julia is the “new kid on the block”: a language which has been designed from the ground up for efficiency and speed, as well as convenience and ease of use. It is supposed to have the power of Matlab, the speed of C, and the elegance of Python. It is still in the early stages (the most recent version is 0.3.0), but it shows enormous promise, and has garnered an impressive list of users and developers. And there is already a small – but growing – list of add-on packages. The Julia web pages show some timings which seem to indicate that Julia is faster (sometimes by several orders of magnitude) than its contenders. (Well, they would, wouldn’t they? They would hardly show timings which were slower than other environments.)

My particular interest is image processing, in which Octave does very well (its image package is quite mature now), followed by Python (with its Mahotas and scikit-image packages). Scilab has its SIP toolbox, but I spent several days trying to install it and gave up. There is an Images package for Julia, but it doesn’t seem on a first glance to have the breadth of the others.

I ran a single test myself: to create a 1000 x 1000 random matrix, invert it, and find the trace of the product of the original and its inverse. In Octave:

tic; A = rand(1000,1000);Ai = inv(A);trace(A*Ai);toc

Scilab:

tic; A = rand(1000,1000);Ai = inv(A);trace(A*Ai);toc

Python (importing numpy as np):

%timeit A = np.matrix(np.random.rand(1000,1000));Ai = A.I;np.trace(A*Ai)

Julia:

tic();A=rand(1000,1000);Ai = inv(A);trace(A*Ai);toc()

Note that Octave, Scilab and Julia all have similar syntax; Python is slightly different because functions don’t all exist in the top namespace. This means, for example, you can use a function called “rand” from different packages, with different functionality, easily.

And here are the results (running the code three times each):

Octave: 1.08677; 1.04735; 1.10616 seconds

Scilab: 3.502; 3.425; 3.663 seconds

Python: 2.33; 2.28; 2.39 seconds

Julia: 0.47034; 0.35776; 0.36403 seconds

This is hardly definitive, but it does show how fast Julia can be, and in relation how slow Scilab can be.

These four environments are by no means the end. R is extensively used, and has over 5000 add-on packages at the CRAN repository. However, I’ve never used R myself. I believe it has branched out from its statistical foundations, and is now more of a general use environment, but still with the best statistical functionality in the business. Then there’s the Perl Data Language, apparently used by astrophysicists, and the GNU Data Language, about which I know nothing. See the Wikipedia page for a list of others.

(Note that I’ve not discussed Freemat this time; I don’t believe it’s a serious contender. I had a look through some of its source code, and most of its files seem to be cut-down Matlab/Octave files without any error checking. Anyway, there seems to be no need to use niche software when Octave is so mature and so well-supported by its developers and users.)

My recommendations? Use Octave, unless you have some niche requirements, in which case you may find Scilab more suitable. If you’re prepared to sacrifice maturity for speed, give Julia a go. If your wish is a mature programming language as your base, then Python.

A very long-running program

In the interests of random number generation, I’ve been experimenting with the iteration

x\to g^x\pmod{p}

for a prime p and a primitive root g\mod p.  Now, it turns out that some primes have a primitive root which generates all non-zero residues, and others don’t.  For example, consider the prime p=19, which has primitive roots 2, 3, 10, 13,14,15.  Starting with x=1, we find the iterations produce:

\displaystyle{\begin{array}{ll}2:&2, 4, 16, 5, 13, 3, 8, 9, 18, 1\\  3:& 3,8,6,7,2,9,18,1\\  10:&10,9,18,1\\  13:&13,15,8,16,9,18,1\\  14:&14,9,18,1\\  15:&15,8, 5, 2, 16, 6, 11, 3, 12, 7, 13, 10, 4, 9, 18, 1  \end{array}}

We see that no primitive root generates all non-zero residues. However, for p=23, we find that the primitive root 20 does generate all non-zero residues:

20,18,2, 9, 5, 10, 8, 6, 16, 13, 14, 4, 12, 3, 19, 17, 7, 21, 15, 11, 22, 1

Let’s call a prime such as 23, an “exponentially complete” prime. You can see a list of the first few such primes at http://oeis.org/A214876 (which has been contributed by me).

My question of the moment is: “Is 2^{31}-1 exponentially complete?” The prime 2^{31}-1 is one beloved by random number theorists, because modulo arithmetic can be performed very efficiently using bit shifts. In C, for example, z \mod p can be computed with

(z&p)+(z>>31)

However, I can’t think of any method of finding this out except by testing each primitive root, to determine the length of its cycle. Even using GNU C, and the PARI library, and some speed-ups (for example, replacing modular exponentiation with just one multiplication from pre-computed arrays), my desktop computer has been grinding away for two days now, and is up to primitive root 1678; that is, the program has taken 2 days to test 380 primitive roots. Since there are \phi(2^{31}-2)=534600000 primitive roots, at this rate it will take something like 534600000/380\times 2 = 2813684.21 days, or about 7703.66 years. That’s a very long time.

The HP Prime Graphics calculator

In my effort to find the perfect CAS calculator, I recently acquired one of these babies:

HP_prime_front_picture

It has many of the hallmarks of the engineering excellence which used to be such a selling point with HP calculators: the keys are superb, the unit itself is thin with a metal back (it has a rechargeable battery, rather than using AAA’s), and its “feel” is of a superior product. The main problem with its externals is the orange letters on the white keys; especially on the numeric keys (which are darker than the others), these letters are hard to read in anything less than very good light. The cyan also on the keys (reached by using Shift) is only marginally better. The keys would be more readable if they were black rather than white or light grey, somewhat like the older HP 35s. HP have always put a lot on to their calculators, in an exemplary fashion, so I don’t know why they missed the mark here.

But that is quibbling – what’s it like for mathematics? In very many ways, absolutely excellent. The screen is crisp, clear, coloured, and also touch sensitive, so that menus can be whisked through by scrolling with the fingers, as well as by using the four-way pad just beneath the screen. My first silly test (factor 2^{67}-1) was done instantly. Very good!

I then decided to see if I could obtain a list of the “zig-zag numbers” z_k; these are the number of permutations of \{1,2,3,\ldots k\} in which the numbers alternately rise and fall. For example, if k=4 these are the zig-zag permutations:

1,3,2,4;\qquad 1,4,2,3;\qquad 2,3,1,4;\qquad 2,4,1,3;\qquad 3,1,4,2

There are five of them, so that z_4=5. It turns out that the exponential generating function for z_k is \tan x+\sec x; so that

\displaystyle{\tan x +\sec x=1+x+z_2\frac{x^2}{2!}+z_3\frac{x^3}{3!}+z_4\frac{x^4}{4!}+\cdots }

You can read about these lovely things on the Wikipedia page on alternating permutations.

So:


s:=series(TAN(x)+SEC(x),x=0,16)
seq(coeff(s,x,k)*k!,k=0..16)

produced the sequence nicely. All of the constructs in those commands are available through menus.

The calculator has two “views”: “Home” and “CAS”. In “Home” view you have floating point calculations, and variables are by default upper-case letters; in “CAS” view operations are exact, and variables are by default in lower-case. This business of lower-case/upper-case confuses me, as I don’t see the need for it. And in fact the two views seem to differ more in their restrictions than anything else. For example, the menus seem the same in both views, but if you try the above series computation in “Home” view, the screen gives you:


CAS.series(TAN(X)+SEC(X),X=0,16)

It looks as though all the CAS operations are available in Home as methods in a CAS class, but in fact the output of this command is simply 1. I don’t understand this at all.

I also use lists a huge amount in my CAS work, but the MAKELIST command only seems to work in Home view, even though it’s available on the same menus in CAS view. Again, I’m quite unclear about the wisdom of having menu items which don’t do anything.

When I introduce students to CAS calculators or to computer-based systems, I get them to find (by trial and error) the smallest k for which

\displaystyle{\sum_{n=1}^k\frac{1}{n}>10.}

This gives them a tiny insight into the speed and power of computational mathematics. And I have done this over the years with Derive, Maple, Maxima, Axiom, Sage, TI-nspire, and CASIO ClassPad. But guess what? you can’t perform this computation on the HP Prime because the SUM function only allows 1000 elements to be summed! I’m told this is an “inadvertent” oversight left over from the HP 39gII, on which much of the HP Prime is based (as well as using a slightly cut-down version of Xcas). It’s a pain though.

The calculator has 14 “Apps”: you can see 12 of them in the picture above (the other two are “Polar” and “Sequence”). In each App some of the keys on the keypad have actions particular to the App. For example, in the “Function” App, the “Symb” key allows you to enter functions, but in the “Sequence” App you get to enter a sequence definition instead. However, the on-screen menus are unchanged.

When I first got the calculator I was plagued by sudden and inexplicable freezes and spontaneous reboots. These have become more rare, but they still occasionally happen. Sometimes the calculator can be unfrozen by a key combination (On+Symb usually works). There is also a “reset” hole at the back, for use with a paper-clip. At the very worst you can unscrew the battery pack cover and temporarily unseat the battery.

I like this machine very much, in spite of its annoyances: I like it a great deal more than the Casio ClassPad (even though the CP has a bigger screen). I would hope that its peculiarities: occasional menu items which don’t work, restrictions on things like SUM, confusions between upper and lower case variables and the Home/CAS views, will be ironed out in the next generation (if there is one).

The Casio Classpad fx-CP400

The Casio ClassPad II (fx-CP400) is the latest CAS offering from Casio: building on the previous ClassPad model to provide a calculator with a large color screen, and the ability to operate it in landscape or portrait mode.

I acquired one of these babies a few months ago, and should have posted about it then, but didn’t. So here goes.

The first reaction on seeing it was that it was very big. In fact I think this is just a matter of its huge screen: the entire unit is in fact not much bigger than the TI-nspire. Here’s a shot of the two of them (in their protective cases), with a pair of spectacles and a ballpoint pen for comparison (the ClassPad is on the right):

cas_calcs1

And here with their cases off, ready for action:

cas_calcs2

The screen is absolutely gorgeous – I can’t say enough good things about it. Crisp, well saturated colors, great for things like sketching a function and its derivatives; really beautiful. What else is good about it? Well, here’s where things get a bit tricky: in fact, I can’t really find any.

Here’s a few random comments:

  1. The system is not much unchanged from the previous ClassPad models. You still have the stylus-driven menus (and very good they are), and woefully underpowered hardware. I did a tiny test: to factorize Cole’s number 2^{67}-1. Here are some timings:
    • TI-nspire: Pretty quick (a few seconds; sometimes longer if the memory’s full)
    • Android Maxima (on my Samsung Galaxy S III): almost instantaneous
    • Casio ClassPad: a day until it ran out of batteries
  2. I have long been critical of 3D graphing on the tiny screens of most CAS calculators: but here at last is a screen on which 3D graphics would work very well. But guess what: Casio, in their wisdom, have not included 3D graphing in this system!
  3. The buttons below the screen are small, plasticky, and feel cheap. Because of their size they wobble a lot. (However, the buttons on the TI-nspire, in spite of comparable size, have hardly any wobble.) The feeling is of a cheap build.

I expect that Casio is aiming fully for the school market. The combination of software and hardware means that this machine would be inadequate for any but the simplest mathematics. I don’t see why this should be so – better programming and a more powerful chip would make this machine into a really useful mathematical tool. As it is I see it as a gorgeous body with not much strength underneath.

A first look a WeBWorK

In my last post (nearly three months ago) I commented on online assessment, and in particular on Pearson’s MyMathLab Global, which we have been using for our engineering mathematics. This is quite a good, full-featured and robust package; its only fault (if indeed a fault it be) is that it’s commercial. This works fine if you are structuring your course around a set text, but over the years I have moved further and further away from such texts. I have not found a text which could be used comfortably with my extraordinarily heterogeneous student cohort, who enter our courses with all different levels of mathematics competency, and all different levels of written and spoken English.

A good text must be carefully, simply, and clearly written, with plenty of useful examples, and above all well scaffolded exercises which move gently from the very simple to the complicated. Most (all?) textbooks in my experience suffer greatly from this. Our current preferred textbook: “Modern Engineering Mathematics” by Glyn James, has quite appalling exercises. Also, it’s bit of a heavy monster of a book, like most mathematics textbooks, and I don’t see why students should have to lug it around with them. But anyway, if you do buy the book, either in hardcopy, or in electronic form, you get a code which gives you one years access to MyMathLab and all ts trimmings.

So this last semester I have used my own notes and exercises for engineering mathematics, and sent students off to MyMathLab for their weekly online tests. There were problems:

  1. Students who first started over a year ago (who, for example, might have had to take a break between their mathematics subjects) found that their access codes had expired. I spent the first few weeks of the semester writing almost daily emails to the publishers begging for some more access codes. To their credit, the local Pearson reps were in every way helpful, accessible, friendly, and understanding.
  2. Right towards the end of the semester I kept finding students who had fallen through the cracks, so to speak, and hadn’t been able to access their tests at all. (And had got busy with other subjects and kept forgetting about them)
  3. Some answers were marked wrong, when they weren’t. I had one student who sent me a screen shot to show that he had entered “12″ as the numerical answer to a problem, which was marked wrong: the correct answer, apparently, was “12/1″. There was also confusion with decimal places: a question would ask for say 4 correct decimal places, and no matter how carefully the students entered their answers they were marked wrong.

I have wanted to experiment with WeBWorK for a while now, but I couldn’t find a way to install it locally – until I realized that I could install it (and a web server) on my desktop at work, which runs Ubuntu 12.04. I did this (it took most of a day, with only a modicum of swearing, and messages to the WeBWorK forums), but now I have a working system, even if only running on a desktop, rather than on a dedicated server.

WeBWorK has many excellent features:

  1. It’s free: open source, no less. You can pay for MAA to host your course if you don’t want to set up your own local service, but even that option is very inexpensive.
  2. It has a huge problem library: the Open Problem Library has something around 25,000 problems from most areas of undergraduate mathematics. You might not find problems on, say, cohomology or modular forms, but these advanced topics would hardly be well served by an online homework system. We are assuming here we have a large class (in the many 100′s) of students in their first or second year of university study. And for this sort of basic mathematics, WeBWorK is terrific.
  3. The authoring system for creating new problems is very straightforward: just a matter of some elementary editing. I tried authoring a problem in MyMathLab once, and although it’s quite possible, life’s too short.
  4. WeBWorK seems to play very nicely with LaTeX and Sage, so that the mathematics is properly typeset, and you can have all the power of Sage available to you if you need it.

That’s as far as I’ve got so far. You can get a feel of WeBWorK without installing it by checking out some of MAA’s “Model Courses”, of which a list is here. There doesn’t seem to be any online hosted method for browsing the problem libraries to see what’s in them, as far as I know. I think this could be something worth making available so that teaching staff can decide whether they think WeBWorK would suit them and their courses.

In my current state of thinking, if I can get a better hosted system at my work (that is, not my desktop), I may well use WeBWorK for my next teaching semester.

Online assessment in mathematics

This semester I’ve taken over a large first year subject, for which the previous subject convenor had organized to use MyMathLab Global for weekly testing.  The subject is based around Glyn James’ “Modern Engineering Mathematics”, a book which is OK for content, and pretty awful for exercises.  This means that as users of the text, we have access to the publisher’s (Pearson) online testing service.  For educators the idea is terrific: every week I simply pull from the exercise bank a set of 10 exercises corresponding to our current topic, and make them available for a week during which time students can have unlimited goes at them.  So in theory it’s an easy way for students to get some easy marks, and reduce the burden of marking weekly tests by the lecturer (me) and the subject tutors.

The subject I am teaching “Engineering Mathematics 2″ is a follow on subject from – wait, you;ve guessed it! – “Engineering Mathematics 1″.

What could be better?  Well, aside from the extraordinary ease of testing, I have begun to have doubts about the efficacy of the system.

  1. It’s a commercial system, which means that students have to buy a “personal access code” to use it.  A code comes with the book if they purchase it.  However, a code lasts for only 12 months, and if students have bought the book (and its code) for Eng Maths 1 (as they all should), and if there’s been a break in their studies for any reason, then their code may have expired by the time they come to Eng Maths 2.  Then there are all the students who have got subject exemptions from Eng Maths 1 and never had a code in the first place.  The local publishers reps have been terrific and have provided me with lots of extra codes to hand out, but the onus is on me to get these dratted codes in the first place, and ensure the students get them.
  2. The system only tests the final answer.  A typical question is:

    Find the partial derivatives

    \displaystyle{\frac{\partial f}{\partial x}}, \displaystyle{\frac{\partial f}{\partial y}} and \displaystyle{\frac{\partial f}{\partial y}} of f(x,y,z)=\sin^{-1}(xyz).

    and the students get three entry boxes for their results.  MyMathLab Global is very picky about the form of a solution, and if the form of a student solution differs from the “correct one” they are marked wrong.  For example,  we find that

    \displaystyle{\frac{\partial f}{\partial x}=\frac{4yz}{\sqrt{1-16x^2y^2z^2}}}.

    if a student leaves out the “1″ of the “16″ in error, the answers is just marked wrong.  If the student decides to write

    \displaystyle{\frac{4yz}{\sqrt{(1-4xyz)(1+4xyz)}}}

    it is marked wrong.  MyMathLab Global doesn’t seem to include a CAS to check if two answers are equivalent (aside from some very simple arithmetic and operator commutativity).  This has been a source of annoyance to my students, who may enter an answer which is mathematically correct, and yet marked wrong.  A very spirited attack on MyMathLab Global can be seen here.

  3. The system checks only the final answer. We might well ask: what mathematical learning is actually being tested here? Surely we want our students to show that they have mastered a mathematical technique, process, or algorithm, and we would mark them on their working, not just on the final answer. At least, that is how I mark by hand. A question for which the student’s working is fundamentally correct but the final result is wrong will still get most marks from me. I regard the final result as just one part of a question. So suppose the student sets out to answer a question, and makes a small slip somewhere along the way which effects every succeeding step, and the final result. I would give almost compete marks – there’s only been one very small slip, and aside from that the logic and the working are exemplary. But MyMathLab Global can’t do this. An answer is either correct, or wrong.

I have been hoping to check out MAA WeBWorK as a comparison, but the only server I have access to runs a linux distribution (RHEL 5.9) which can’t run it. I think WeBWorK, from my understanding, is rather more nuanced in terms of mathematical equivalence of expressions, so probably overcomes item 2 above, and being open source completely overcomes item 1. I don’t think it overcomes item 3, though.

I will continue to use MyMathLab Global – its convenience is simply too great – in spite of my misgivings. But I think that online assessment in mathematics is still a very long way from testing the students’ ability to do mathematics, as opposed to their ability to obtain an answer. This latter is not mathematics: it is applied arithmetic; and noble and useful as it may be, it is only a tiny part of what we, as mathematics educators, should be testing.