Calling the WMAP likelihood code from C/C++

If you’re interested in cosmological parameters, chances are you’ll want to include WMAP CMB constraints in your parameter estimation code at some point. Happily, the WMAP team have made their likelihood code publicly-available, so it’s relatively simple to add them into your own MCMC software (or whatever else you’re using). Less happily, the likelihood code is written in Fortran 90, so users of more modern languages will need to do a bit of fiddling to get it to play ball with their code.

Joe Zuntz was kind enough to share a little C wrapper that he wrote for the WMAP likelihood code. It’s pretty easy to follow, although you should of course read the documentation for the likelihood code to see how to use it properly. First of all, compile the original Fortran WMAP likelihood code. Then, compile this wrapper function as a static library (libwmapwrapper) as follows:

gfortran -O2 -c WMAP_likelihood_wrapper.F90
ar rc libwmapwrapper.a WMAP_likelihood_wrapper.o

You can call the wrapper function from your C code by linking to that library. Joe has written a test implementation that shows how it works. To compile this, you’ll need to make sure you’re linking against everything the WMAP code needs, including a BLAS/LAPACK library; the following should work on Mac OS X (using veclib):

gcc wmap_test.c -std=c99 -L. -lwmap -lwmapwrapper -lcfitsio -framework veclib -lgfortran -o test_wmap

(N.B. Joe’s code is written for the WMAP 7-year release, so you may need to change a couple of numbers to get it working with the 9-year release.)

About these ads

About Phil Bull

I'm a postdoctoral fellow in cosmology at the Institute of Theoretical Astrophysics, Oslo. My research focuses on the effects of inhomogeneities on the evolution of the Universe and how we measure it. I'm also keen on stochastic processes, scientific computing, the philosophy of science, and open source stuff. View all posts by Phil Bull

3 responses to “Calling the WMAP likelihood code from C/C++

  • Phillip Helbig

    “Less happily, the likelihood code is written in Fortran 90, so users of more modern languages will need to do a bit of fiddling to get it to play ball with their code.”

    Why the problem with Fortran? The newest standard is Fortran 2008 (what is the latest C or C++ standard?), which actually includes a standardized C interface, but I seriously doubt that Fortran 90 is not sufficient for the task at hand. It is a subset of more modern Fortran. So, why cast Fortran as old-school? Much code, even in cosmology, is written in it.

    Similarly, many people no longer work on VMS, where it is common practice to write an application in a variety of languages. (Indeed, VMS is written in a variety of languages.) C and C-like languages are today’s monocultures. Be afraid of the next storm.

    • Phil Bull

      My problem with Fortran is that it makes it easier than most languages to write nasty, messy, hard-to-maintain code (Perl has similar issues; it’s often described as a “write-only” language). That’s not to say that all Fortran code is nasty and messy — in the right hands, it can work well. But unfortunately, past experience shows me that Fortran doesn’t find its way into the right hands that often…

      Some languages do a lot to promote good programming practices (Python being a good example). While you can write crappy code in any language, some make it more difficult. Compared with Python or C++, I’ve found Fortran 90 often seems to make it easier to do things the “wrong way” (e.g. it’s easier to use global variables and goto statements) than the “right way”. Unless you’re a Fortran guru, or have beaten the bad habits that Fortran allows out of yourself, the result is often spaghetti code — spaghetti code that I sometimes end up adopting, and thus resenting.

      Fortran 90 is about as old-school as C in terms of syntax and language features (array and string handling is easier, but that’s about it). Fortran 2008 has more modern features — I’m thinking of things like OOP — but it’s not so widely used or supported. I don’t want to say that newer necessarily means better, but there are languages out there with some really wonderful new(ish) features and syntactic conventions that save me a lot of time (both in writing and debugging), Python being my favourite. Of course Fortran shouldn’t be written off — it has its place — but in general, I think other languages can do what Fortran does in more succinct, and less error-prone, ways, making them better choices for new projects.

      (P.S. The affinity for C-like languages is probably something to do with the dominance of UNIX-based systems, which are very C-based.)

      • Phillip Helbig

        Thanks for the balanced reply. You have a point. :-) I think this is due mainly to Fortran being the first high-level language and code lasting decades so that often some legacy stuff is still around which could be done better with a rewrite (but is still legal for backwards compatibility). As someone once said, Fortran programs aren’t written, they are rewritten.

        Someone also said that a Fortran programmer can write Fortran in any language. :-) However, the important point is that one can write good code in modern Fortran, either by knowing what one is doing (always a good thing when programming) or by using some subset of Fortran which has been developed for this purpose.

        Nevertheless, I often find bad Fortran more legible than good C. :-(

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: