Matlab vs. Python

I've never programmed a single line of Python code but recently a lot of people have been telling me that it has a lot of neat advantages over Matlab. I think I'm going to have to find a day to learn Python's ins and outs and get a feel for when it makes sense to use it over Matlab. This article was quite persuasive to me. Time for a look. I'll report back once I've tried out some basics.

Numerical Simulation

The other day I came across a neat problem that I'd used to teach the basics of numerical simulation to some undergraduates a few summers ago. I thought I'd re-post it here. The goal is to solve for the solution to the following third order differential equation:'

• y''' + 5y'' + 8y' + 4y = 1

Use a timestep of dt = 200us and run the simulation for t = 10 seconds. The initial conditions are:

• y(0) = 5
• y'(0) = 1
• y''(0) = 0

This problem is fun because you can solve for a closed form solution and then compare that hand-solved answer to the simulated version. The simulated solution can be run using Matlab, C, or a combination. The combination solution is neat because it uses Matlab to create the data and plot the solution, but uses C-code (compiled in Matlab as a mex-function) as the super-fast solution engine. Mex-Functions are a bit tricky to learn but can often lead to valuable speed-ups in simulations.

The solution to the third order differential equation can be solved by hand. You should get the following function:

The plot for this function is shown below along with the iterated solution that I generated using C++:

Simulated Solutions

Here I present three numerical (i.e. simulated) solutions to the third order differential equation. All three yield the same simulated plot (see blue curve above).

Matlab-Only Solution

clear; 

clf;

dt = 200e-6;

tMax = 10;

nSamples = tMax/dt;

y = zeros(nSamples,1);

a = 1;

b = 0;

y(1) = 5;

for i=2:nSamples

    dy = a;

    da = b;

    db = -8*a - 5*b - 4*y(i-1) + 1;

    y(i) = y(i-1) + dy*dt;

    a    = a      + da*dt;

    b    = b      + db*dt;

end

t = (1:length(y))*dt;   

plot(t,y);

Using the tic/toc commands in Matlab, I determined that the Matlab-only solution took an average of 12ms, including the time to create the plot.

C++ Only Solution

The C++ Only solution is fast, but we need a way to pass the data back to Matlab in order to plot it. In this solution, C++ writes the data to a binary file. Then a separate Matlab script reads the data from the file and plots it.

#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;
int main()
{
double dt = 200e-6;
        double tmax = 10;
double pi = 4 * atan(1);
int    nSamples = floor(tmax/dt);
double da,db,dy;
double a=1;
double b=0;
double y[nSamples];
int i;
ofstream out("data2.bin",ios::out|ios::binary);
y[0] = 5;
for(i=1;i<nSamples;i++){
dy = a;
da = b;
db = -8*a - 5*b - 4*y[i-1] + 1;
y[i] = y[i-1] + dy*dt;
a    = a      + da*dt;
b    = b      + db*dt;
}
out.write((char *)&dt , sizeof(double));
out.write((char *)y  ,nSamples*sizeof(double));
out.close();
return 0;
}

This is the Matlab plot code:

clear;clf
fid = fopen('data2.bin','rb');
dt  = fread(fid,1,'double');y   = fread(fid,'double');
fclose(fid);
t = (0:length(y)-1)*dt;
plot(t,y);xlabel('time (s)');title('Solution to 3rd Order Diff-Eq');

C++ / Matlab / Mex Solution

While the previous solution works, it requires that we switch back and forth between the Matlab and C++ environments, which is inherently inefficient; it also requires that we generate a large data file for the purpose of shuttling the data back and forth. A better solution is to create a mex-file. A mex file is a C/C++ file that is compiled directly within Matlab. The compiled executable can be called directly from Matlab; data parameters can be passed back and forth from Matlab to the executable without the intermediate step of dumping it in a file. This solution requires that we only work with one programming environment: Matlab. The coding is a bit more complicated, but the solution is ultimately more elegant. The complexity of the coding is primarily due to the way Matlab passes data to C/C++. The Matlab data comes in structures, with pointers everywhere; learning to maneuver in this manner takes some getting used to. There are good references for this process here and here.

#include <mex.h>
#include <string.h>
#include <math.h>
// This is the subroutine that actually performs the simulation
void runSim(double **py, double **pt, double dt, double tMax, int *nSamples){
    double pi = 4 * atan(1);
    double da, db, dy;
    double a=1;
    double b=0;
    int i;
    double tTemp = dt;
    double *y, *t;
   
    *nSamples = floor(tMax/dt);
    *py = new double[*nSamples];
    y = *py;
   
    *pt = new double[*nSamples];
    t = *pt;
   
    y[0] = 5;
    t[0] = dt;
    for(i=1;i<*nSamples;i++){
        dy = a;
        da = b;
        db = -8*a - 5*b - 4* y[i-1] + 1;
       
        y[i] = y[i-1] + dy*dt;
        a    = a      + da*dt;
        b    = b      + db*dt;
       
        t[i] = t[i-1] + dt;
    }
}
// ****************************************************
// ******************** START HERE ********************
// ****************************************************
// Mex routines must always start with "mexFunction"
// Here, the input data is imported from Matlab. Then the actual function
// executed (in this case "runSim"), and finally the output data is
// exported back to Matlab
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    // Step 1: Import input variables from Matlab
    double dt   = *(double *)mxGetData(prhs[0]);
    double tMax = *(double *)mxGetData(prhs[1]);
    // Step 2: Declare variables
    int nSamples;
    double *y,*t;
   
    // Step 3: Run the simulation
    runSim(&y, &t, dt, tMax, &nSamples);
   
    // Step 4: Export the results back to Matlab
    double *output;
    if (nlhs>=1){
        plhs[0] = mxCreateDoubleMatrix(nSamples, 1, mxREAL);
        output = mxGetPr(plhs[0]);
        memcpy(output, y, nSamples*sizeof(double));
    }
    if (nlhs>=2){
        plhs[1] = mxCreateDoubleMatrix(nSamples, 1, mxREAL);
        output = mxGetPr(plhs[1]);
        memcpy(output, t, nSamples*sizeof(double));
    }
   
    // Step 5: Housekeeping
    delete [] y;
    delete [] t;
} 

Once the C/C++ code has been written and compiled, the function is called from Matlab in the same way that any other function would be. In this case, the C file was called "diff_eq3.cpp", so therefore the compiled executable is called using "diff_eq3". This simulation took an average of 2ms, including plotting time. That is 6x faster than the Matlab-only solution!

clear; clf;
dt = 200e-6;
tmax = 10;
[y,t] = diff_eq3(dt,tmax);
plot(t,y); 

Dennis Ritchie

Unix and C pioneer Dennis Ritchie has passed away. Bummer. This IEEE Spectrum article gives the full details on his amazing career. I'll remember him for being one of the authors of the world's best C programming book, best known as "K+R" or "Kernighan and Ritchie". I had the pleasure of hearing Kernighan speak at Temple University a couple of years ago... pretty cool stuff. Anyways, thanks for the phenomenal tool, Dr. Ritchie. You've probably touched more lives that even Steve Jobs did!

Formatting Cout

I just found a neat solution to a common but vexing problem: formatting cout in C++. Turns out its pretty simple. This sample program shows how to enforce two decimal places of precision:

---

#include <iostream>
#include <iomanip>

using namespace std;

int main()
{
double x = 800000.0/81.0;
double y = 1.2345;

cout << setiosflags(ios::fixed);
cout << setprecision(2);

cout << x << endl;
cout << y << endl;

return 0;
}

---

>>
9876.54
1.23

---

The first formatting command enforces a fixed-point output, and the second enforces two decimal places of precision.

A much more complete explanation of options can be found at this blog.

Neat Coding Trick # 3

Here's the final useful coding trick I learned last week. While Matlab is unparalleled for rapid prototyping, there are any number of reasons why you might want to deploy your code in another language such as C/C++. For example, you might be targeting an embedded device that requires C/C++. Or you might just want to port and compile certain functions you've written in Maltab. This would be the case if you've written a function that is slow to execute in Matlab but that will run very quickly once compiled. Just about anything that is loop-intensive will fall into this category.

For a while now I've known about "mex" functions. Mex functions are C/C++ routines which you compile at the Matlab command line; mex functions must include a special function called mexFunction which specifies how data is passed back and forth between Matlab and the C routine. Once the code is compiled, you call the function from Matlab just as if it were a Matlab function and hopefully you realize huge savings in execution time. This technique works pretty well but the downsides are (a) You have to write the code yourself from scratch and (b) you have to write the mexFunction interface which is tricky at best.

Happily, I discovered "emlc"which generates C code from Matlab code. Here's how it works.

1. Write a function in Matlab
2. Use emlc to convert that code to C from the command line.

Wasn't that easy? There are a few important options that are worth being aware of. Suppose your function is "myFunc.m". You could go to the command line and type:

• emlc -T MEX myFunc

In this case,  your myFunc.m would be converted into myFunc.mexmaci64 (or some other extension depending on your operating system). Your code is automatically ported over and compiled. So you can do all your prototyping in Matlab and then use this simple command to get a high speed executable.

If you are targeting an embedded system and you just need the raw C code, you can instead try using:

• emlc -T RTW

There are a host of options; I recommend reading the emlc documentation to be aware of all the options.

Here are a couple of useful tricks I've learned by trial and error. Most importantly, whereas Matlab will dynamically allocate memory for you on the fly, C/C++ is not set up to do the same and must be explicitly instructed on how to handle memory. For example, suppose you've written a function that takes as input an array "y" which will have two elements in it. Whereas Matlab will give you the benefit of the doubt that "y" will in fact have two elements at run time, C will not. To avoid a compile-time error, you must include the following statement just after your Matlab function heading:

• assert(all(size(y)==[2,1]));

This command tells the compiler to throw a run-time error if the "y" input to the function is not exactly of dimensions two rows and one column.

Along the same lines, you must also pre-allocate all output variables. If your output variable is "dy" and you expect it to have two elements, you can't just say:

• dy(1) =  ... ;
• dy(2) =  ... ;

Instead, you have to say:

• dy = zeros(2,1); % preallocate
• dy(1) =  ... ;
• dy(2) =  ... ;

I also learned that its not good enough to preallocate with an indefinite number. For example, this will cause a compile-time error:

• function dy = myFunc(y,tmax)
• int nPts = tmax/0.001;
• dy = zeros(nPts,1);

The problem is that the compiler has no idea how many points nPts will be. The workaround is to preallocate up to some maximum worst-case array length and then to truncate it as necessary:

• function dy = myFunc(y,tmax)
• int nPts = tmax/0.001;
• int ptsMax = 1e6;
• dy = zeros(ptsMax,1);
• dy(1) = ... ;
• dy(2) = ... ;
• dy = dy(1:nPts);

Finally, you need to include the following comment next to your function heading: %#eml. This facilitates the creation of the compilation report:

• function dy = myFunc(y,tmax) %#eml

Just for kicks, I coded up a single-cell Hodgkin-Huxley simulation using Forward Euler. The point was to create something slow and lumbering so that I could really see the speed difference. You can download the mFile by clicking here.

I used the following command at the Matlab prompt:

• clear; tic; [v,t] = hodgkin_huxley(10,1,100); toc; plot(t,v);

Over five trials, the average elapsed time is 62ms.

Then I converted the m file to a mex file using the steps outlined above and reran it using exactly the same command line command:

• emlc -T MEX hodgkin_huxley
• clear; tic; [v,t] = hodgkin_huxley(10,1,100); toc; plot(t,v);

Now the average elapsed time is 14ms. Thats an execution-time reduction of over 75%, and all it cost me was typing in a single line a the command prompt.

Neat Coding Trick # 2

My Day of Coding last Friday turned up a really useful trick for passing data between Matlab and C/C++ during execution. Suppose you have some C++ code and lots of data, and you want to visualize or otherwise examine your data during runtime either to debug your code or just to make sure that things are happening the way you think they are. Matlab has provided tools to allow you to do just that: pass data from C/C++, during runtime, to Matlab. There, the data can be visualized or otherwise manipulated, and even sent back to C++. Not too bad.

The basic steps are:

• include engine.h - this is an include file that automatically installs with Matlab
• define an "engine" pointer, which opens up a portal between C++ and Matlab
• define Matlab-friendly data elements (mxArrays) for porting your data to Matlab
• stuff data into the Matlab-friendly data elements
• instruct C++ to send the Matlab-friendly data elements to Matlab
• instruct C++ what commands to execute on the data

There are of course more sophisticated options; see Matlab online documentation for more details.

I've put together a brief video outlining the technique. The code shown in the video can be downloaded from here: http://www.box.net/shared/xmazi14uv2

Neat Coding Trick # 1

Last Friday, I learned three neat coding tricks which I'm going to try to share here. The first one involves inherited classes in C++ and the use of virtual functions to create run-time decisions about which versions of a function to run. This is a very powerful trick because it simplifies function calls and data handling in situations where you have a number of related classes.

I've put together a little demonstration here to illustrate how useful this technique is:

---

#include<iostream> 
using namespace std;

class base {
protected:
 int x;
public:
 void setX(int val){x = val;};
 int getX(){return x;};
 virtual void incrX() = 0;
};

class derived1 : public base {
public:
 void incrX(){x+=1;};
};

class derived2 : public base {
public:
 void incrX(){x+=2;};
};

int main(){

 base *myVars[2];

 myVars[0] = new derived1;
 myVars[1] = new derived2;

 myVars[0]->setX(0);
 myVars[1]->setX(0);

 myVars[0]->incrX();
 myVars[1]->incrX();

 cout << myVars[0]->getX() << endl;
 cout << myVars[1]->getX() << endl;

 delete myVars[0];
 delete myVars[1];

 return 0;
}

---

We can think of this code as implementing two versions of a common class. The "common" portions of the class are in the base class. The base stipulates a private variable named x, as well as functions for setting and retrieving x. The base also stipulates an undefined virtual function named incrX. The way I've coded, this, any derived class that inherits base must define its own implementation of incrX. In the "derived1" class, incrX increments the value of "x" by 1, whereas in "derived2", incrX increments the value of x by 2.

In the main function, I create an array of pointers of type *base. I can then create new instances of the derived class and have those pointers stored in the array of type *base. This is a pretty amazing trick. Because derived1 and derived2 both inherit base, I am allowed to define a pointer of type *base and point it to either of the derived classes.

The second remarkable part of this code is that when I run the myVars[0]->incrX() command, the code is clever enough to realize that myVars[0] actually points to an object of class "derived1"; it then runs the appropriate version of incrX.

This trick is very handy because it is going to allow us to solve a few nasty problems in our neuron simulator. We have models for about five basic neuron types. In many ways, those neuron types are similar: all must keep track of who their pre-synaptic neurons are and all must have a function for numerically updating the state variables. However in other ways, those neurons are quite different: the differential equations and state variables are all different from neuron to neuron.

The elegant solution to this model is to create a "base" neuron which contains all the elements that are common to all neurons. The base neuron will also stipulate a virtual "update" function which will need to be redefined by each specific neuron type. Then we can create five "derived" classes which inherit the base and add the individual update functions and state variables. The great part is that in the "main" function, I only have to maintain a single array of neurons. I do this by creating an array of pointers of type *base. Then I can populate that array with any combination of the five neurons. When I tell a neuron to "update", the program makes sure that the update function appropriate to the specific neuron is called.

Another neat upshot of this technique is that the base neuron class can contain a vector of pointers to base which can be populated with pointers to presynaptic neurons, regardless of their specific type.

Thanks to Chris for figuring most of this out!