Kirua Studio

Math Love

2009-12-02T18:31:00.003+01:00

Matlab: Tips from the Trenches

2009-09-13T20:05:00.010+02:00

Every once in a while, you'll stumble upon something you can't do in Matlab, and spend way too many hours figuring it out. Every once in a while, you will spend hours doing things manually because you didn't even suspect Matlab could do it for you. Here are a few tips from the trenches that could save you a lot of time and help you build better software.

Often times, you will be writing a LaTeX report to communicate your Matlab simulation results. This gets tiresome if every time you modify your data or simulation parameters you have to regenerate and save every plot and array by hand to replace it accordingly in your LaTeX document. We'll see how all of this can be automated by learning to use a few commands. Remember to type 'help ' at the Matlab command prompt to find out more.

Perfect plots all by script, and how to get them in LaTeX

The first command you really should know about is 'print'. It lets you save the current figure to an image file from script or the command prompt. I have two favourite modes:

print('-dpng', '-r100', 'filename.png');

This will save the current figure in a PNG image file named 'filename.png', with a resolution of 100 dots per inch. This is useful for quick communication or backup of the plot, since PNG's are widely supported. My other favourite mode is:

print('-depsc', 'filename.eps');

This will save the current figure to a colour EPS file. This is awesome because EPS images are vector based, hence: whatever the resolution you print them, they will look perfect. And LaTeX can exploit that feature! Basically, an EPS file contains all the data Matlab used to draw the plot in the first place, and enough instructions to redo the job for any display size or paper printing. Unfortunately, EPS files are less convenient for communicating your results as such, since only few people know how to display them and have the proper software handy.

Usually, you will want to use both commands, to keep a quickly accessible version of your plot, and a neat version for your LaTeX report. So your code would look something like this:

plot( ...... );
print('-dpng', '-r200', 'myplot1.png');
print('-depsc', 'myplot1.eps');

Let's find out how to include the EPS graphics in a LaTeX document that we will export to PDF (works running Windows for sure, should be just the same running Linux). First, you need to include the graphicx package:

\usepackage{graphicx}

Then just include your figure like you would always do, only make sure to *not* specify the extension of the image file name. So for an image called filename.eps, write something like this:

\begin{figure}[h]
\centering
\includegraphics[width=.75\textwidth]{filename}
\caption{Description of my perfectly drawn plot.}
\label{myfigurelabelforreferences}
\end{figure}

The trick now is to *not* call PDFLatex. The latter cannot deal with EPS files. You need to first build your LaTeX source using the normal LaTeX compiler, which will output a DVI. Then call DviPDF on that file. An example of a fairly good free LaTeX editor running on Windows that lets you do that easily is LEd (LaTeX Editor). That's it!

Now, when you save plots automatically, you want to control every single detail of how it will appear from script, since you won't be able to do it manually using the GUI. It's actually something you want to do anyway, since the GUI can get very slow. Here are a few tips on how to do just that (call matlab help on them):

figure : will open a new window for plotting purposes

title, xlabel, ylabel, zlabel: work as usual to set the title and axis labels. You might want to set the font size, example:

xlabel('\fontsize{12} Frequency \xi [Hz]');

Indeed, you can use greek letters with the usual LaTeX commands.

xlim, ylim, zlim : choose the limits of the axis box

Here is how you can choose where the axis ticks are.

set(gca, 'XTick', 0:10:100);

Here, 'gca' means "get current axes", which lets you work on the current figure, really. XTick means you want to act on the X axis, but you could use YTick or ZTick. The last parameter is a vector indicating where you want the axis ticks to be.

grid on or grid off displays or hides a grid over the plot, according to the axis ticks.

This lets you set the axis ticks font size too:

set(gca, 'fontsize', 12);

Here is a very handy function, that lets you impose the aspect ratio of the different axes:

pbaspect([1.3 1 1]);

This means that the X axis will be 30% longer than the Y axis (and the Z axis is not important for a 2D plot). You can also use:

axis equal : means that a unit on the X axis will span the same size (in centimeters) as a unit on the Y axis. That comes in handy if the X and Y axis have the same units.

axis tight : prevents the plot box from having a white unused border region: the box will be defined by the extreme values of your plotted data.

One more thing: if you happen to work in a display less environment (we'll talk about this some more later on), you can still use all of these functions, and that's a pretty convenient way of still being able to produce plots. You can also force Matlab to not display a figure, so that generating bunches of plots in a loop would be faster:

set(gcf, 'visible', 'off');

Here, 'gcf' means "get current figure", and you should call this function after you created the figure by calling the 'figure' command.

Porting Matlab code to C++

From time to time, you'll need to do this. It's a pain. Really. But less so if you use the right tools. I gave the GSL a try: do NOT use that library. Complex numbers and matrices support (or lack thereof) are awful. These are the tools I use and enjoy:

Armadillo For vectors and matrices support: very high quality code, with delayed evaluations: incredibly efficient, amazing syntax, motivated developer. Highly recommend it. It integrates with LAPACK and BLAS (and ATLAS and BOOST). If you're a windows user, you might as well not lose to much time trying to compile these libraries (it's horrible, take my word for it), and get these binaries from Victor Liu. You will need a latest generation GCC compiler for Armadillo to work. The neat thing is: Matlab is built on BLAS and LAPACK, so you will get the same numerical results.

FFTW is the FFT computation library of choice, and is also the one that Matlab uses. It integrates just fine with Armadillo, with a little care. This library is pretty sophisticated, so you might want to spend some time reading the documentation.

For the rest, you're on your own. I'm not going to lie to you: it will be painful.

Matlab via SSH and the screen command

If you get access to a computation server via SSH (that's pretty common), you might want to run long simulations on the server, without the need to keep the SSH connection alive. Here is how you can do that:

'ssh' into your server
execute the 'screen' command (just type 'screen' and hit return)
launch matlab in console mode: 'matlab -nodisplay'
hit 'ctrl+a' then 'd' (this will detach you from the screen)

You're now back to your previous terminal, and matlab disappeared. How can you get back to it? Just type 'screen -r' from *any* terminal (on *any* computer!) that is connected to your server session via ssh. This will get you back to visualizing what your Matlab program is doing.

This way, you can launch a computation via ssh from your laptop, then shut your laptop down, go back home, get some sleep, come back the next day and just get your results back. If you want to monitor the activity of your process, the 'top' command might be useful (inside or outside the screen program, it doesn't change anything). Screen is a very powerful command, and there is much more to it than just this.

A word about transpose and hermitian

I was shocked to discover that, after four years of intense Matlab use: the " ' " (quote) operator is *not* transpose! It's conjugate-transpose (hermitian). The proper operator for the transpose is " .' " (dot-quote). I didn't know about that important difference, and found out that many of my friends didn't either. I dearly hope it was just us.

Profiling your Matlab code and memory management

Sure, Matlab is only for prototyping. But most of the time, it's just not worth your time to port a Matlab code to C or C++, and you still need to run pretty big computations. So how can you speed things up? Well, my advice would be to use the profiler to spot the tricky parts of the code, and write some C-MEX code for these bits. The profiler basically monitors how much time each line of your code has been executed, which lets you find out what bits of code make up for most of the execution time. And so you know where to focus your attention to get things going faster.

Here is basically how to use it:

profile clear;
profile on;

% execute your code here

profile off;
profile viewer;

A GUI will pop up and give you plenty of information about your code.

Remember though: optimization is the last step in software development. Early optimization is pretty much an excuse for code obfuscation. A readable, documented, maintainable, straightforward, simple code is what every one wants. When you have that, profile. The profiler will tell you where it's worth it to do some dirty weird coding.

Also, please note that the 'reshape' command is actually free: no memory is copied around or moved on the physical medium. That is because all matrices and arrays are stored as 1D vectors. Reshaping a matrix really just changes the way you access elements with multi indices.

An important rule to remember about matlab memory management is this:

When you call a function, every parameter passed to it is, by default, passed as a reference. Only if the parameter is modified inside the function will it be passed by copy.

This means that it *is* okay to pass a 1 Gb matrix as a parameter to a function that you call a billion times. As long as the data in the matrix is just read, and not modified, it will not be duplicated in memory: no overhead.

Example:

function s = getsize(A)
s = size(A);
end

A will not be copied: free.

function s = getsize(A)
A = A + 1;
s = size(A);
end

If you have doubts about whether your data is being passed by reference or by copy, here is a nice way to find out. Type this at the command prompt:

format debug;

Now, every time you display a variable, Matlab will also display the memory pointer to the variables data. It's called the 'pr' field. Here is an example:

>> A = rand(3);
>> B = A;
>> A

A =

Structure address = 8fb42e0
m = 3
n = 3
pr = c2d8a90
pi = 0
0.7060 0.2181 0.3709
0.6451 0.7724 0.8909
0.5523 0.2280 0.8564

>> B

B =

Structure address = 8fb7c60
m = 3
n = 3
pr = c2d8a90
pi = 0
0.7060 0.2181 0.3709
0.6451 0.7724 0.8909
0.5523 0.2280 0.8564

As you can see, A and B share the same 'pr' field: A and B point to the *same* memory. Now if we actually modify B, Matlab will copy A's data and save it in a separate location:

>> B = B + 1

B =

Structure address = 8fb7c60
m = 3
n = 3
pr = adbb2e0
pi = 0
1.7060 1.2181 1.3709
1.6451 1.7724 1.8909
1.5523 1.2280 1.8564

Use the same principle with function inputs / outputs to make sure that your data is not being stupidly copied around. But Matlab does a pretty good job, so not to worry. To revert to normal display again, type 'format short'.

That's it for today,

Have fun optimizing the world!

Kirua

Trivial

2009-07-28T11:36:00.003+02:00

United Against

2008-12-27T14:05:00.001+01:00

Harmonic Inversion

2008-12-26T18:16:00.003+01:00

Just an update for the last post, even though it's really old now:

This library, HarmInv, from MIT, would have done the frequency estimation job just right. Too bad I didn't know about it. It actually decomposes a signal into a finite sum of exponentially decaying sinusoids, without doing a least mean squares regression (which could lack robustness, since it's not a linear problem in this case). Brilliant idea!

Since this summer, I had some image processing classes, and I can think of a million other ways of doing this over, harder, better, faster, stronger. I'll probably come back to it later. Specifically, the Hough transform was a nice idea to play around with, but it's really overkill and some form of prediction (using a Kalman filter for example) would have been nice.

Measuring Stuff with a Webcam

2008-09-13T18:01:00.014+02:00

Long time no write!

I set up an original experience in my bedroom. It's designed to measure a length using a webcam. On this picture, you can see how things are laid out:

The principle is quite simple:

build a pendulum the same length as the object you want to measure,
point a webcam in front of it,
connect the webcam to your computer,
set the pendulum in motion, and acquire a bunch of frames,
find the pendulum position in each frame,
knowing the physics of the situation, analyze the motion to determine the length of the pendulum.

It's easy because:

you only need to evaluate the oscillation frequency: it's no big deal if the camera is not set up perfectly parallel to the motion plane or if the resolution is poor,
mass, initial conditions and amplitude (if tiny) are nothing to worry about,
the webcam delay is not a problem, since only the time difference between two frames matters,
you don't need to convert from cartesian coordinates to polar ones, since you know the actual motion is close to a decaying sinusoid: counting zero-crossings is sufficient.

It's not so easy because:

Using a simple model for the pendulum (massless rod, infinitely small bob), you don't compute the correct length. If we choose to neglect moment of inertia though, which I did, the length you measure is the distance between the pivot point and the center of mass of the pendulum. Hence, you should choose a fairly small bob, which I couldn't (because I'm too clumsy and can't fix a string on a tiny object).

I used Matlab and the excellent Image Acquisition Toolbox to acquire frames from my cheap webcam. Here is a snapshot:

The next step is to locate the pendulum in this image. I wrote a modified Hough Transform to identify circles in order to do that. To convert the disc to a circle, I used a gradient filter, on the green channel to emphasize the edges:

Once and for all, I manually estimated the radius of the pendulum. Running an HT on this yields the following Hough-space (ie: parameter space):

Extracting the maximum intensity in the above image gives the cartesian coordinates of the pendulum position. I did this on a sequence of 3000 frames, captured at 25 FPS with a resolution of 320x240 pixels. I pinched the pendulum beforehand, not caring whether it was in the right alignment etc: it doesn't matter. Here is an excerpt of the data, for the horizontal position (vertical position doesn't vary enough to be relevant):

As expected, the observed "signal" is a decaying sinusoid. We wish to evaluate the frequency of the latter. IMHO, a Fourier analysis isn't the best choice here, since we know that the signal much ressembles a stable carrier amplitude modulated by a very low frequency decaying function. Simply counting zero-crossings is acceptable in this configuration. According to this formula:

l = g/w^2

where l is the pendulum length, g is 9.81 m/s^2 and w = 2*pi*f, where f is the oscillation frequency, the test yields the following results:

f = .958 Hz
l = 27.1 cm

It's funny that I unwillingly built a nearly 1 Hz pendulum ^_^.

Of course, the measured length covers the distance between the pivot point and about the center of the bob.

Currently, the Hough Transform is too slow, so the computations are done offline. To speed things up, one could use a model of the system to estimate the next position of the bob, and only look it up in a small area.

The initial project was to mount the pivot point on a horizontal axis. The position of the pivot would then be controlled by a motor, itself controlled by the Lego Mindstorms RCX brick. In order to do so, I made the necessary physics calculations and coded an RCX Matlab Toolbox (interesting to do actually), but the latter was way too slow to expect anything interesting to be feasible. Damn the infrared protocol ... The idea was to pinch the pendulum, and have the computer automatically move the pivot point in order to stabilize it as fast as possible. It's a nice project, because it involves a lot of math, a bit of physics and practical issues. Unfortunately, it's not realizable with the RCX. At least, not with my software.

If you're aware of any way to control the motors of the RCX brick live from a computer, with a reasonable delay, in whatever language, I would be glad to hear about it!

Kirua

XKCD is Bad for You

2008-07-10T20:20:00.002+02:00

You might want to take a look at this: xkcd nerd sniping.

I confronted a handful of my math and electricity teachers with this problem during classes. It's funny how they usually freeze for a couple of minutes, unable to resist the temptation to try and solve it. My linear control teacher even spent a full 15 minutes not giving class to think about it. My favourite math teacher told me to come back the next week, which I did. He simply walked to me (ooh, how I hate it that I was sleeping during his class at that very moment), and quickly said: it has an elegant solution, then left. He wouldn't tell me more. I ended up solving the problem empirically with a Matlab simulation but ... it simply doesn't taste like victory at all.

Fundamental Laws

2008-07-10T18:42:00.003+02:00

Long time no see! There should be more on the way.

Fives Stages

2008-06-08T12:40:00.004+02:00

Aside from the beer, the fifth stage also involved improvising lots of songs on the theme "internet is down" and being mocked by our roommate.

Automatic Caller Detection Circuit

2008-05-10T16:37:00.003+02:00

For some couples, you might need specially designed small capacitors to build a very high pass filter. In some extreme cases, a normal microphone won't do. I recommend getting those Coustaud used to record whales.

Food for Thoughts

2008-05-10T15:26:00.006+02:00

Whatever.

Speaking about rockstars and food for thoughts, here are two quotes I'd like to share with you (my keyboard?) :

Life's a journey, not a destination.
Aerosmith

Someday, somehow, I'm gonna make it alright, but not right now.
Nickelback

I Think We're Ready

2008-05-03T18:36:00.005+02:00

Partial credits to Karim

If I Had Written Die Hard

2008-04-07T22:20:00.004+02:00

Beta-hydroxybutyryl-CoA dehydrogenase. It hurts.

The Pendulum Equation

2008-03-28T23:59:00.007+01:00

True Random Numbers

2008-03-26T23:33:00.003+01:00

These great kids at the INMA ...

By the way, if you're going to look that acronym up, make sure to enter "INMA UCL" and *not* just "INMA".

Last week, we got our new applied math sweat-shirt. It's red and there's an xkcd comic in the back. It bears a graph-like logo, and we figured only on the delivery day that it had an eulerian path. How cool is that!?

...

Yeah, exactly!

Untitled

2008-03-24T13:16:00.003+01:00

Cartoon GUI

2008-03-23T22:56:00.003+01:00

Yayks! I built myself a Matlab GUI for my photo to cartoon conversion algorithm :). This way, I can test multiple tuning parameters very easily and cut out parts I don't need in the picture. Love it :). The last two posts were made using this program.

Pancakes

2008-03-23T22:40:00.005+01:00

Looking at this plot, it seems impossible to conceive a rationnal chain of events that could have lead to the invention of pancakes. Everybody in his right mind would have stopped at the batter stage. The inventor of pancakes did a major "outside the box" thinking.

Today, Mr. Pancake is chairman at the optimization and operations research department of my university. It's a well earned honour sir, way to go!

Credits to Karim for this one.

Speed Dating

2008-03-23T22:20:00.003+01:00

Apple Crumble Recipe

2008-03-14T22:35:00.003+01:00

It's hard to tell the complexity. But if execution time is an issue, you can very simply parallelize this algorithm (in case you've got a multi-core cooking team). More refined versions implement a negative feedback loop inherited from control theory in order to regulate the appropriate baking time, but that would need a good quality sensor, like a girl for instance, who would be capable of telling whether the crumble needs more baking or not. The experimental constant used here should provide sufficiently accurate results though.

Caution! Make sure not to include this routine in a loop! Stomachs have been known to segfault.

Next time, I'll expose a couple of optimization algorithms to improve the input quantities and other parameters, resulting in optimal satisfaction.

P=NP Lab

2008-03-11T00:07:00.005+01:00

Just in case, I also drew this one for a proof by contradiction and an other one using the pizza corollary. You just never know, it might come up. Lately, my friends and I have gotten quite confident that our professors do know that P=NP and just won't tell. Every once in a while, we set out to try and confund them. We'll keep you tuned.

By the way, that's another one drawn on a blackboard. It's really nice actually, way easier than on regular paper. Too bad I didn't plan the content enough, though. The "storyline" doesn't feel quite right.

Applied Automatics

2008-02-25T22:27:00.006+01:00

G's output variable is the volume of fluid in the bottle. The input variable is the angle of the bottle, with respect to the horizontal, defined by the gravity vector g. The idea is as follows: every person in a party is an autonomous system which has it's own drinking dynamics, it's own behaviour. By putting people in a closed control loop though, you can force them to drink much faster, effectively building "servo-drinkers".

Science. It works, bitches! ... And it makes great parties, too ^^.

Just for the record, the cartoon was actually drawn on a blackboard with an actual chalk. It's kind of a tribute to my algorithm that nothing needed to be changed in the script I usually use to convert the following picture to the above cartoon:

Plenty of ideas for the next ones, so stay tuned!

All Apologies

2008-01-31T18:26:00.001+01:00

That's it for today! Six cartoons in a row after more than a month drawing nothing. I've got plenty of ideas, some starring Matlab as a hero, but these are yet to be refined.
See you soon, hopefully ^_^.

Malevich

2008-01-31T18:24:00.000+01:00

Theology Today

2008-01-31T18:17:00.001+01:00

Asking his students about uncomputable problems, here is what a teacher from my university's computer science department got as an answer:

determining the last decimal digit of pi
dividing by zero
counting to infinity ... twice

Including Chuck Norris in today's computability theory definitely is this decade's challenge.

Shrupman's Filter

2008-01-31T18:09:00.001+01:00

I'm getting really interested in stochastic control theory these days. I got this book for christmas, and it's a delight. Yet, it would be unreadable without a course about stochastic processes, which I had the chance to follow last semester. Unfortunately, the Automatics course I choose for this semester seems to be somewhat loosy... Wait and see.

Godwin's Law

2008-01-31T18:06:00.001+01:00

(This is a redrawing of an old one)

Happy New Year 2008!

2008-01-31T18:02:00.001+01:00

This trick will stop working counting from 2010 for a long time. Next year's is a fine one too, though:

(2+0!)!*(0!+sqrt(9)!)

Charging your iPod with Onions?

2007-11-19T23:39:00.000+01:00

Yesterday evening, I saw this video on the net:
"How to Charge an iPod using electrolytes and an onion"...
http://www.youtube.com/watch?v=GfPJeDssBOM

I got obsessed by the idea. There's something really strange in this procedure: no electrodes! For an oxydo-reduction to take place in the onion, there should be electrodes, right?

So I collected a multimeter, a big onion, lots of salt, some water and a screwdriver (most of which was readily available in my student appartment), and I followed the instructions (except for the last iPod part, and I used salted water as an electrolytic solution instead of Gatorade or Powerade which were suggested).

Well I'm telling you ... I've got my doubts!

Voltages across the onion were really small (they quickly decreased from about 70mV to less than 10mV). I gathered two coins (made of two different but unknown metals) and plugged them into the onion, then mesured voltages across them: voltages went as high as 250mV! And longer, too. It's been five hours now, and the voltmeter still displays 35mV. It was still at more than a hundred only 90 minutes ago.

Now, as far as I know, USB connectors are made of a single metal (probably copper or gold, I wouldn't know). And I don't believe oxydo-reductions take place on them anyway ^^.

So ... What's the trick here?

I'm certainly no chemist, so if I'm missing something, please tell me. Because I'd certainly want this video not to be a fake! Thinking of the impression it might make just drives me mad ^^.

Timing Distances for a Living

2007-11-19T23:11:00.001+01:00

This one is a memory from my trip to Dubrovnik, Croatia. The meeting of two maths engineers, a physicist and a geographer. Surprising combination. 'specially for the geographer one might say ^^.

I've been timing heights since I'm 15 or so. We went from timing our distance to the sea on the city walls to using that data to determine the distance to the horizon.

Quite funny, really:

h = gt²/2
d = sqrt(2Rh) (assuming h << R)

So (approximately):

d = sqrt(Rg) t = 7.91 . t

where d is expressed in km and t in seconds.

Too bad air friction plays such an important role here, even when there is little or no wind. Results are horrible ^^.

The next step involved timing height by jumping in the sea. But these pictures you won't (and don't care to) see!

Natural-Behaviour Optimization Algorithm

2007-11-18T23:43:00.000+01:00

Simulation results were not so bad actually.

I'm not too sure what impact this post could have on future job interviews ...

Euclid 1, Kolmogorov 0

2007-11-11T01:14:00.001+01:00

The comic is supposed to be on a single line, but I had to rearrange it for this blog.
I initially planned an alternative ending where the statistician wins... Then I changed my mind.

Is math invented or discovered?

2007-11-11T00:26:00.000+01:00

If Mathematics were invented and not discovered,
then why did Earth orbit around the sun
before Mathematicians were born?
In a universe without Mathematicians,
does a tree falling in an empty forest make noise?

From the moment you define the axioms of a theory, every single theorem that holds in it... exists. It doesn't wait for you to imagine it. Accepting that, I have some trouble with the idea that math could have been invented. Then again, maybe this question is more about the concept of axiom and theorem rather than about the actual ones.

Something is rotten in Chemistry

2007-11-09T19:53:00.002+01:00

What else can I say ... Chemists won't even admit there's a problem here!

My Pick-up Line

2007-11-09T19:46:00.001+01:00

Well, of course, it does narrow your chances of getting laid. You're making a really good point here...

Belgian Political Issues

2007-11-09T19:39:00.001+01:00

To people outside Belgium: please bear in mind that nobody's fighting for real here. These are all political issues which do not prevent people from both sides to actually get along together great. Well, most of them though.

Oh, and, of course: you can assume that the bridge is massless for further calculations.

Random Joke

2007-11-09T19:36:00.000+01:00

Lorem ipsum dolor sit amet, consectetuer adipiscing elit. Nunc ut urna quis tellus interdum venenatis.

Murphy's Law Restated

2007-11-08T23:17:00.001+01:00

This one opens a series, I guess. It's not the first one I drew, but it applies so well to reality. By the way, how do you define reality?

Reality is what doesn't cease to exist when you stop believing in it.
Philip K. Dick

Open bar, counting from now.

2007-11-08T22:52:00.000+01:00

Hello potential virtual readers!

This is the official first post of this blog, which will be mainly about my cartoons, usually drawn during long boring speeches of the very few less gifted teachers at my university. Which basically means there will not be many of them.

The intent is to be funny ... that you can certainly contest and comment, I'm sure. Due to my lack of any talent in that area though, the aesthetic aspect of cartooning will not be the point here (sic).

The geek thing about all this really is ... that the algorithm I use to convert my poor quality noisy pictures to black and white png's was developped at night with matlab (god I love that stuff!), based on mainly original ideas. The final tuning was inspired by my all time partner though: wikipedia.

Special credits go to my geeky community: Coder-Studio.com.

I hope you'll enjoy reading the sometimes nonsensical stories here!

Kirua, a.k.a. Nicolas Boumal.

PS: by the way, this material is very clearly inspired by the great XKCD webcomic, but by no means tries to compare to it.