The Chromatic Fantasy and Fugue in D minor, BWV 903 by Johann Sebastian Bach played by Canadian pianist Angela Hewitt, perhaps the best Bach interpreter since the great Glenn Gould.

## Selection of the week

September 19, 2014## Julia update

September 16, 2014I’m totally committed to Julia now. It is blitzing fast and very similar to Matlab with some small differences (improvements). I particularly like the fact that when you declare an array with one argument, like x = zeros(10), it immediately gives you a vector of length 10 and not a 10 X 10 matrix. The broadcast function is also very useful. There are still many things I don’t yet know how to do in Julia like how to import and export data. Plotting is also not fully solved. I’ve been using PyCall to import matplotlib.pyplot, which works pretty well but not perfectly. There are things I miss in Matlab like context dependent history recall, i.e. I can recall an old command line by just typing the first few letters. If anyone knows how to do this in Julia please let me know. Right now, I’m hitting the uparrow button continuously until I find the line I want. I do worry that my abandoning Matlab means someone will lose their job, which is certainly not my intention. However, I am well ahead of schedule for zero Matlab.

## Selection of the week

September 12, 2014I think the Meditation from Jules Massenet’s opera Thais is appropriate for the day after September 11. Here is violinist Sarah Chang.

## Unintended consequences

September 10, 2014Here is a true story. A young man is trained to hit people as hard as possible and to react immediately to any provocation with unhindered aggression. He signs a 5 year contract for 35 million dollars to do this 16 times a year or more if he and his colleagues are very successful at doing it. One day he gets upset with his fiancée and strikes her in the head so hard that she is knocked unconscious in a public place. This creates a minor stir so the employer mandates that he must apologize and is prohibited from smashing into people for 2 of the 16 times he is scheduled to do so. The fiancée-now-spouse also refuses to press charges because she doesn’t want to jeopardize the 27 million over the next 3 years owed to the man. However, a video showing the incident is made public creating a huge uproar so the employer abruptly fires the man and condemns him since he now is no longer financially useful to the employer. The public now feels vindicated that such a despicable man is no longer employed and that domestic violence now is given the attention it deserves. However, the spouse is very unhappy because her comfortable lifestyle has just been pulled from right under her. Now, other spouses who live with violent but rich men will be even more silent about abuse because they fear losing their livelihoods too. If we really cared about victims of domestic violence, we would force the employer to set up a fund to ensure that spouses that come forward are compensated financially. We would also force them to support institutions that help the many more victims of domestic abuse who are not married to rich and famous people. This young man is probably an upstanding citizen most of the time. Now he is unemployed and potentially even angrier. He should not be thrown out onto the street but given a chance to redeem himself. The employers and the system who trained and groomed these young men need to look at themselves.

## The morality of watching (American) football

September 7, 2014The question in this week’s New York Times Ethicist column is whether it is wrong to watch football because of the inherent dangers to the players. The ethicist, Chuck Klosterman, says that it is ethical to watch football because the players made the decision to play freely with full knowledge of the risks. Although I think Klosterman has a valid point and I do not judge anyone who enjoys football, I have personally decided to forgo watching it. I simply could no longer stomach watching player after player going down with serious injuries each week. In Klosterman’s article, he goes on to say that even if football were the only livelihood the players had, we should still watch football so that they could have a livelihood. This is where I disagree. Aside from the fact that we shouldn’t have a society where the only chance to have a decent livelihood is through sports, football need not be that sport. If football did not exist, some other sport, including a modified safer football, would take its place. Soccer is the most popular sport in the rest of the world. Football exists in its current form because the fans support it. If that support moved to another sport, the players would move too.

## Selection of the week

September 5, 2014I think we should start September off with Beethoven. Here is the first movement of the Sextet in E-flat op.81 b for two horns and strings performed by the Amici Ensemble Frankfurt. Note that they added a double bass, just to make it more special.

## Selection of the Week

August 29, 2014Chopin Nocturne No 8 Op 27 No 2 performed by Maurizio Pollini, whom I’ve had the fortune of seeing in Boston many years ago.

## MCMC for linear models

August 25, 2014I’ve been asked in a comment to give a sample of pseudo code for an MCMC algorithm to fit a linear model to some data. See here for the original post on MCMC. With a linear model, you can write down the answer in closed form (see here), so it is a good model to test your algorithm and code. Here it is in pseudo-Julia code:

# initial guess for parameters a and b

a=0

b=0

# construct chi squared, where D is the data vector and x is the vector of the

# independent quantity

chi = norm(D - (a*x +b))^2;

for n = 1 : total;

# Make random guesses for new normally distributed a and b with mean old a and b

# and standard deviation asig and bsig

at = a + asig * randn()

bt = b + bsig * randn() chit = norm(D - (at*x + bt))^2;

# Take ratio of likelihoods, sigma is the data uncertainty

ratio=exp((-chit + chi)/(2*sigma^2));

# Compare the ratio to a uniform random number between 0 and 1,

# keep new parameters if ratio exceeds random number

if rand() < ratio a = at;

b = bt; chi = chit; end

end

# Keep running until convergence

## Selection of the week

August 22, 2014George Gerschwin, although generally classified as jazz, was strongly influenced by classical music. Here is Summertime from his opera Porgy and Bess performed by soprano Renee Fleming.

## J. Bryce McLeod, 1929-2014

August 20, 2014I was given the sad news that J. Bryce McLeod died today in his home in England. Bryce was an extraordinary mathematician and an even better human being. I had the fortune of being his colleague in the math department at the University of Pittsburgh. I will always remember how gracious and welcoming he was when I started. One of the highlights of my career was being invited to a conference in his honour in Oxford in 2001. At the conference dinner, Bryce gave the most perfectly constructed speech I have ever heard. It was just like the way he did mathematics – elegantly and sublimely.

## Julia vs Python update

August 19, 2014We’ve been using Julia for a little over a month now and I am quite pleased thus far. It is fast enough for us to get decent results for our MCMC runs; runs for which Python and Brian were too slow. Of course, we probably could have tried to optimize our other codes but Julia installs right out of the box and is very easy to program in. We still have issues for plotting but do have minimal plot functionality using PyCall and importing matplotlib.pyplot. One can also dump the results to a file and have either Python or some other software make the plots. I already find myself reaching for Julia when I want to write a quick code to solve a problem instead of relying on Matlab like I used to. While I would like to try PyDSTool, the lack of a trivial installation is holding us back. For you budding software developers out there, if it takes more than one click on a link to make it work on my computer, I am likely to go to something else. The reason I switched from Linux to Mac a decade ago was that I wanted Unix capability *and* the ability to print without having to consult with a sys admin.

## Selection of the week

August 15, 2014Eiene Kleine Nachtmusik by Wolfgang Amadeus Mozart performed by the Concertgebouw Kamerorkest in 2013 in Amsterdam.

## Numerical integration

August 14, 2014Solving differential equations numerically is a very mature field with multiple algorithms carefully explained in multiple textbooks. However, when it comes down to actually solving them in practice by nonspecialists, a lot of people, myself included, will often resort to the Euler method. There are two main reasons for this. The first is that it is trivial to remember and code and the second is that Moore’s law has increased the speed of computing sufficiently that we can get away with it. However, I think it’s time to get out of the eighteenth century and at least move into the nineteenth.

Suppose you have a differential equation , where can be a scalar or a vector in an arbitrary number of dimensions. The Euler method simply discretizes the time derivative via

so that at time step is given by

Now what could be simpler than that? There are two problems with this algorithm. The first is that it is only accurate to order so you need to have very small steps to simulate accurately. The second is that the Euler method is unstable if the step size is too large. The way around these limitations is to use very small time steps, which is computationally expensive.

The second most popular algorithm is probably (4th order) Runge-Kutta, which can generally use larger time steps for the same accuracy. However, a higher order algorithm is not necessarily faster because what really slows down a numerical scheme is the number of function evaluations it needs to make and 4th order Runge-Kutta needs to make a lot of function evaluations. Runge-Kutta is also plagued by the same stability problems as Euler and that can also limit the maximum allowed step size.

Now, these problems have been known for decades if not centuries and work-arounds do exist. In fact, the most sophisticated solvers such as CVODE tend to use adaptive multi-step methods, like Adams-Bashforth (AB) or it’s implicit cousin Adams-Moulton and this is what we should all be implementing in our own codes. Below I will give simple pseudo-code to show how a 2nd order Adams-Bashforth scheme is not much more complicated to code than Euler and uses the same number of function evaluations. In short, with not much effort, you can speed your code immensely by simply switching to AB.

The cleverness of AB is that it uses previous steps to refine the estimate for the next step. In this way, you get a higher order scheme without having to re-evaluate the right hand side of your ODE. You simply reuse your computations. It’s the ultimate green scheme. For the same ODE example above 2nd order AB is simply

You just need to store the previous two evaluations to compute the next step. To start the algorithm, you can just use an Euler algorithm to get the first step.

Here is the 2nd order AB algorithm written in pseudo Julia code for a one dimensional version of the ODE above.

# make x a two dimensional vector for the two time steps

# Set initial condition

x[1]=0;

# Store evaluation of function f(x,t) for later use

fstore[1]= f(x[1],h)

# Take Euler step

x[2] = x[1] + h*fstore[1]

# store function at time = 2h

fstore[2] = f(x[2],2*h)

# Store computed values for output

xout[1]=x[1]

xout[2]=x[2]

# Take AB steps

for t = 3:Tfinal

# set update indices for x based on parity of t

index2 = t%2+1

index1 = (t-1)%2+1

# 2nd order AB step

x[index1] = x[index2]+h*(1.5*fstore[2]-0.5*fstore[1]))

# update stored function

fstore[index1] = f(x[index1],t*h)

xout[t]=x[index1]

end

return xout

## Selection of the Week

August 8, 2014The 4th movement of Franz Schubert’s Trout Quintet (Piano Quintet in A Major) performed by the Schubert Ensemble.

## The neuro software gap

August 4, 2014Our foray into Julia was mostly motivated by a desire to have a programming language as fast as C but as easy to program as Matlab. In particular, we need really fast looping speed, which is necessary to simulate dynamical systems. Our particular project involves fitting a conductance-based model (i.e. Hodgkin-Huxley-like equations) to some data using MCMC. We’re not trying to fit the precise spike shapes *per se* but the firing rates of synaptically coupled pools of neurons in response to various inputs. We thus need to run simulations for some time, compare the averaged rates to previous runs, and iterate. For good fits, we may need to do this a million times. Hence, we need speed and we can’t really parallelize because we need information from one iteration before going to the next.

Now, one would think that simulating neural networks is a long solved problem and I thought that too. Hence, I instructed my fellow Wally Xie to look into all the available neural simulators and see what fits the bill. We first tried Neuron but couldn’t get all of it to work on our Linux machine running CentOS (which is what NIH forces us to use). That was the end of Neuron for us. We next tried Nest but found it difficult to specify the neuron model. We finally settled on Brian, which is adaptable, written in Python, and supported by very nice people, but is turning out to be too slow. Native Python is also too slow. Stan would be great but it does not yet support differential equations. Wally is now implementing the code in Julia, which is turning out to be very fast but has a possible bug that is preventing long runs. The developers are extremely helpful though so we believe this will be worked out soon. However, it is hard to get plotting to work on Julia and as pointed out by many, Julia is not nearly as complete as Python or Matlab. (Yes, Matlab is too slow). We could use C and will go to it next if Julia does not work out. Thus, even though people have been developing neural simulators for decades, there is still an unmet need for something adaptable, fast, easy to use, and runs out of the box. A neuro-Stan would be great.

## Selection of the Week

August 1, 2014Summer from Antonio Vivaldi’s Four Seasons of course! Here is a performance by the 2nd place winner at the 2012 Menuhin Junior violin competition.

## JSMB/SMB talk

July 29, 2014I just gave an invited plenary talk at the joint annual meeting of the Japanese Society of Mathematical Biology and Society of Mathematical Biology in Osaka Japan. I talked about my work on steroid-regulated gene transcription. The slides are here. Previous posts on the topic, including background summaries, can be found here.

## Tim’s Vermeer

July 28, 2014Jan Vermeer has been one of my favourite painters ever since I saw his famous “The Girl with the Pearl Earring” painting that was on display in Toronto in the 1980’s. I’ve been on a quest to see all of his paintings although its been on hiatus for the past ten years. Here is the list of what I’ve seen so far (I have five left). You only need to stand in front of a Vermeer for a few seconds to be mesmerized. I stood in front of “The Music Lesson” in Buckingham Palace for at least an hour. The guard started joking with me because I was so transfixed. This is why I’ve been intrigued by recent suggestions by artist David Hockney and others that some great old masters like Vermeer and van Eyck may have used optical aids like camera obscura. Well, inventor Tim Jenison has taken this theory to another level by attempting to completely recreate Vermeer’s Music Lesson using a set up of mirrors and lenses that he (re)invented. The endeavor is documented in the film Tim’s Vermeer directed by Teller of Penn and Teller fame. Whether you believe the theory or not (I actually do and it doesn’t detract at all for my love of Vermeer), what this film does do so well is to show what dedication, thought, patience, and careful execution can accomplish. I got tired just watching him paint the threads in a Persian rug using his optical tool.

## Selection of the week

July 23, 2014What is called classical music mostly refers to the Western symphony orchestra tradition that starts in the seventeenth century with Vivaldi and peaks in the early twentieth with Mahler. While classical music remains popular, my unscientific sampling of concert hall audiences indicates that the demographic skews to retirement age and above. I don’t know if this means that a generation of music lovers is about to depart or that people only have the patience to sit through a long concert when they are older. In an attempt to introduce a new generation to classical music, I thought I would present a selection each week. And what’s a better way to kick it off then with the pseudo-Baroque precursor to heavy metal, the Praeludium and Allegro by Fritz Kreisler. Kreisler performed in the first half of the twentieth century. He was one of the greatest violin virtuosos of all time and also wrote some great violin ditties. Here is a performance by the then 13 year old Canadian/American violinist Leila Josefowicz in 1991.

## Talk at Jackfest

July 16, 2014I’m currently in Banff, Alberta for a Festschrift for Jack Cowan (webpage here). Jack is one of the founders of theoretical neuroscience and has infused many important ideas into the field. The Wilson-Cowan equations that he and Hugh Wilson developed in the early seventies form a foundation for both modeling neural systems and machine learning. My talk will summarize my work on deriving “generalized Wilson-Cowan equations” that include both neural activity and correlations. The slides can be found here. References and a summary of the work can be found here. All videos of the talks can be found here.

Addendum: 17:44. Some typos in the talk were fixed.

Addendum: 18:25. I just realized I said something silly in my talk. The Legendre transform is an involution because the transform of the transform is the inverse. I said something completely inane instead.