Category Archives: Research

Peer-to-peer distributed Matlab computing

In a recent meeting with the SPM developers, we discussed parallel computing using the Matlab distributed computing toolbox, Star-P, Sun Grid Engine, and other batch systems that can be linked to Matlab. These are all limited in their usefulness for the typical neuroimaging research setting in that they are based on a centralized job distribution system. That may work fine on a large cluster with a centralized configuration and system administration, but even then the usefullness is limited because all input and output data (which are typically large) have to be send over the network twice: first to the job manager, then to the compute node (and vice versa for the results).

To resolve some of these problems, I came up with the idea of peer-to-peer distributed computing in Matlab. The full description can be found on

Automatically compile a missing mex file on the fly

I am maintaining a few open source Matlab toolboxes. Although most functions are plain Matlab, some functions are implemented as mex file. The mex files have to be compiled on each platform, i.e. 32 and 64 Linux, 32 and 64 bit Windows, ppc and intel Mac OS X, etc. Since I don’t have access to all those platforms, I have devised a trick to automatically compile the mex files on the end-user’s platform.

Let’s say that you have a function mymexfile, then the file mymexfile.c would contain the c-code implementation of the Matlab mex file. You should save the following code as mymexfile.m. If the compiled mex file (e.g. mymexfile.mexw32) does not exist upon the first call to the function, the m-file will be executed. The m-file will automatically compile the c-file and will subsequently call the compiled mex file. On all subsequent calls, the compiled mex file will be executed directly.

function [varargout] = autocompile(varargin)

% AUTOCOMPILE compile the missing mex file on the fly

% remember the original working directory
pwdir = pwd;

% determine the name and full path of this function
funname = mfilename('fullpath');
mexsrc = [funname '.c'];
[mexdir, mexname] = fileparts(funname);

% try to compile the mex file on the fly
warning('trying to compile MEX file from %s', mexsrc);
success = true;

% compilation failed
error('could not locate MEX file for %s', mexname);
success = false;

if success
% execute the mex file that was juist created
funname = mfilename;
funhandle = str2func(funname);
[varargout{1:nargout}] = funhandle(varargin{:});

Multithreading in a MATLAB mex file

Together with colleagues from the F.C. Donders Centre and the NICI, I am currently working on a project on real-time analysis of EEG and MEG data for Brain-Computer Interfacing (BCI). Since we already have a lot of relevant code in MATLAB and since a lot of the target end-users are familiar with it, we would like to do the real-time processing in MATLAB.

The idea is to read the data that is streaming from the EEG amplifier, do the computation on the data (e.g. spatial filtering, spectral decomposition and classification) and then continue with the next section of the data. The challenge is not to lose any of the streaming EEG data while MATLAB is busy with the computation. A potential solution would be to implement the reading/acquisition of the data into matlab in a mex file which is multi-threaded.

For this purpose I wrote some code for testing the multithreading in MATLAB. The test code below demonstrates that it works, and shows how easy it actually is.

I am working on an Apple PowerBook G4 with OS X 10.4.11 and with MATLAB Version 7.4 (R2007a). The example code uses the POSIX Pthread libraries and I would expect the same code to work on most Unix and Linux systems. It will not work on a windows system, but might work with pthreads-win32.

Put the following code in a file with the name threads.c

#include < pthread.h > /* for threading */
#include < unistd.h > /* for sleep */
#include "mex.h"
#include "matrix.h"

void *do_thread(void *threadid)
int tid;
tid = (int)threadid;
mexPrintf("In thread: just started thread #%dn", tid);
sleep(4); /* pretend that the code is busy for a few seconds */
mexPrintf("In thread: about to exit from thread #%dn", tid);

void mexFunction (int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[])
int num;

if (nrhs<1)
mexErrMsgTxt("not enough input arguments");
num = mxGetScalar(prhs[0]);

pthread_t threads[num];
int rc, t;
for(t=0; t < num ; t++){ mexPrintf("In main: creating thread %dn", t); rc = pthread_create(&threads[t], NULL, do_thread, (void *)t); if (rc){ mexErrMsgTxt("problem with return code from pthread_create()"); } sleep(1); /* wait some time before making the next thread */ } return; }

Subsequently, in MATLAB you should compile the mex file using the following command:

>> mex threads.c -lpthread

After that, you can run the mex function. It takes a single argument, corresponding to the number of threads to create. Each thread will live for three seconds (pretending that it is busy for some time), and the time between threads being created is one second. Some example output is below.

>> threads(6)
In main: creating thread 0
In thread: just started thread #0
In main: creating thread 1
In thread: just started thread #1
In main: creating thread 2
In thread: just started thread #2
In main: creating thread 3
In thread: just started thread #3
In thread: about to exit from thread #0
In main: creating thread 4
In thread: just started thread #4
In thread: about to exit from thread #1
In main: creating thread 5
In thread: just started thread #5
In thread: about to exit from thread #2
>> In thread: about to exit from thread #3
In thread: about to exit from thread #4
In thread: about to exit from thread #5

You can see each thread starting and stopping. The last threads stop after the control has been returned to the MATLAB prompt.

Here you can find some comments from Mathworks, including some windows specific example code. Important to note is that MATLAB itself is single-threaded, and not thread safe. Hence you should take care in interfacing your threaded c-code in the mex file with the MATLAB core.

I would appreciate it to get some feedback on this. Please leave a comment below if you tested it successfully or if it does not work for you, with some details on your configuration (like operating system and MATLAB version).

Average reference for dipole fitting

Here is a question that I get asked occasionally. I have slightly edited the question and my answer to it, which both were posted to the EEGLAB email discussion list.

Question: I’m using EEGLAB-DIPFIT to localize independent components using the spherical head model. Apparently the software requires the data to use the average reference. Why is this?

Answer: In principle you could use an arbitrary reference in your source reconstruction. The practical reason to use an average reference over the sampled electrodes in source estimation is that this prevents the solution to be biassed due to forward modelling errors at the reference electrode. Let me give a partially intuitive, partially mathematical explanation.

Assume that you would use left mastoid as reference. That would mean that the measured value “V” at each electrode “x” is V_x, so the list of all measured values in the N channels is

V_M1-V_M1 (this is zero)

Those values can be modeled using the source model and the volume condution model. Now, lets assume a spherical volume conduction model. That is especially inaccurate for low electrodes, and the bony structure of the mastoid is definitely not modelled appropriately in a spherical model. So for the model potential “P” we would have the value at each of the N electrode also referenced to the model mastoid

P_M1-P_M1 (this is zero)

The source estimation algorithm tries to minimize the quadratic error between model potential distribution and the measurement, so the error term to be minimized is
= sum of quadratic error over all channels
= [(V_C3-V_M1)-(P_C3-P_M1)]^2 + ….
= [(V_C3-P_C3)-(V_M1-P_M1)]^2 + …. (here the terms are re-ordered)

So for each channel the error term consists of a part that corresponds to the potential at the electrode of interest, plus a part that corresponds to the reference electrode. The error term corresponding to the reference electrode is identical over all channels (i.e. repeats in each channel), hence for each channel you are adding some error term for the reference electrode. Therefore, the minimum error (“minimum norm”) solution will be one that especially tries to minimize the model error at the reference electrode (since that is included N times). In the case of a mastoid reference we know that there is a large volume conductor model error at M1, hence the source solution would mainly try to minimize that error term. The result would be that the source solution would be biassed, because it tries to reduce the (systematic) error at the reference.

The solution is to use an average reference (average over all measured electrodes). That implicitely assumes that the model error over all electrodes is on average zero, hence the minimum norm solution is not biassed towards a specific reference electrode.

PS the maths in my explanation above are rather sloppy, but the argument still holds for a more elaborate mathematical derivation which would assume the forward model inaccuracies are uncorrelated over electrode sites.

High-density EEG electrode placement

Some time ago we wrote an article on electrode placement for high-resolution EEG measurement (referred to as the 5% article). After its apearance I have noticed that there is a demand for a concise and methodological overview of electrode placement systems. With this page I want to share some of my knowledge on this subject. This page contains non-technical comments on the different standards for electrode placement. Continue reading

EEGLAB workshop

From 17-19 September 2005, the second EEGLAB workshop will be held in Porto, Portugal. I will be giving a lecture on the use of the DIPFIT plugin (i.e. dipole analysis of ICA components). The new version 2 of the DIPFIT plugin will be introduced during he workshop, which means that the link between Fieldtrip and EEGLAB is finally official. For more information about the workshop, you can look here. For more information about the relation between FieldTrip and EEGLAB, you can look here and here.

Matlab compiler

Currently I am testing out the use of my mentat toolbox for parallell computing within Matlab in combination with GridEngine on our linux cluster. Initially it worked fine, but that was because I was testing it with simple level Matlab functions like “rand” and “fft”. However, I am now trying to get it to work with more complex real-world scripts for EEG and MEG analysis that use the FieldTrip toolbox and I am running into a few problems:

I had to edit the file /opt/matlab-6.5.1/toolbox/matlab/iofun/private/readavi.m and rename the output argument “varargin” into “varargout”. I think that this is a plain bug in that function, which only becomes apparent after translating the code to c. The original file resulted in a compilation error, changing the output variable name fixed that.

Somehow, the code wants to include the Matlab function wk1read in the stand-alone executable. I cannot determine where that happens, because all functions are very deeply nested: mcc generates ~350 c-files on which my analysis function depends. The problem is that, after compilation, wk1read results in a unresolved symbol error. Since I am not interested in reading Lotus database stuff (which is what that particular function does), I have simply replaced it on my own Matlab search path with an alternative (empty) implementation for wk1read. That solved the unresolved symbol error.

Furthermore, I encountered a problems with the exist function. In the stand-alone executable, it does not seem to detect the existence of a variable. Since I am using that function on multiple locations throughout the code, I have to think of a nice workaround. Related to this problem I found this page which gives some relevant hints on using the Matlab compiler:

  1. The following things cause a warning in compilation and may cause a runtime error if the offending statement is executed: objects, most toolboxes, inline, whos, some graphics functions that are passed functions by string.
  2. Program compiled with graphics will not run without X windows. This includes running on beowulf nodes from PBS scripts.
  3. Program containing graphics commands and compiled without graphics will produce compile time warning but will run fine as long as the graphics commands are not executed.
  4. The following commands work differently: eval, feval, exist, global.
  5. Use handles to pass functions rather than strings.
  6. Do not use exist to check for the presence of arguments in a function.