Who needs five digits of accuracy anyway

At the ANS Mathematics & Computational topical conference this spring, I heard a rumor (from a very well-placed source) about an ironic error in MCNP.

MCNP is considered by many in the labs and industry to be "the" nuclear simulation program of choice. It's a Monte Carlo transport program with continuous energy data and the ability to represent virtually exact geometry (with no truncation error for smooth surface, unlike in deterministic methods). It's had many decades of development, and it's been tested and tweaked considerably to match up with experimental and numerical benchmarks. Naturally, then, a simulated result from MCNP is often used as a reference answer when comparing new codes.

No simulation can exactly represent reality: there are always errors in data, simplifications in the physics, etc. This should always be kept in mind when doing computational anything. Case in point, I heard of a "run-off" benchmark between MCNP and another reactor code, SCALE. The result was that MCNP got the right answer and SCALE was wrong: the only problem was that their definition of the "right" answer was based on an MCNP-run calculation.

The error in MCNP—which you can verify yourself with the card entry print 98—is laughably simple yet indeterminately old. Avogadro's number, which converts between the number of atoms and the mass in a material, is programmed to be 6.022043447E+23 rather than the correct 6.02214179E+23. (In fact, it's written to 16 digits of "accuracy" even though only about ten are scientifically known to any precision.) The irony is that if they change Avogadro's number, their carefully benchmarked answers will come out "wrong."

It's not a huge problem, since I doubt most answers would be sensitive to such a small perturbation, but it serves as a warning to carefully consider the difference between "true" and "best approximation."

Grammar fact of the day

I am apparently overly prone to hyphenating words. I knew that compound adjectives are supposed to be hyphenated. (Compound adjectives are modifiers composed of two words stuck together.) However, I wasn't sure if predicate adjectives were supposed to be hyphenated. I would certainly have a time-dependent simulation, but could I say that my simulation is time-dependent? As it turns out, I cannot grammatically say so. My problem is time dependent, as is my comprehension of the English language and its quaint usage rules.

SWIG and Traits classes

With all the new code I've been writing—namely, a large number of classes and template parameters and typedefs—I've been using a number of traits classes. This unfortunately posed a problem for SWIG (which provides a Python front-end to the C++ code), because it had trouble deducing the true object types through the maze (albeit short) of inheritance, typedefs, and template parameters. It had trouble despite the %include "Traits.hpp" command, which makes it aware of the traits class with template parameters.

My first solution was to wrap the troublesome traits-derived typenames with #if SWIG [explicitly substitute types] #else [use normal traits] #endif. It nullified a subset of the purpose of the traits class, but it did the job well enough.

The problem cropped up again today with a new class. Frustrated, I cast about for a solution, and tried explicitly instantiating the template in SWIG with %template(TraitsT) Traits<OneD, Gray>;. To my surprise, it worked! So apparently it's necessary to instantiate a traits class (even if it is a struct with no methods or data) in order to get SWIG to properly parse the typedefs in classes that use those traits.

EDIT: Looks like the manual actually has an entry about this. It recommends %template() traits<double,double>; to avoid generating unnecessary wrappers for it.

Bad writing

I am compelled to link to How To Write Badly Well, whose example texts read like entries to the Bulwer-Lytton Fiction Contest. Coincidentally, they both resemble Dan Brown's writing style.

Make dylib paths absolute

One problem I have come across in my trek through the murky swamps of open source libraries and so forth have been occasional issues with shared libraries on the mac. Shared libraries (dylib files on the mac, so files on Linux systems) are modules of code loaded at run time; the binary executables (or other object files) only store a path to the shared library file (and sometimes a version number). Sometimes this path is wrong: either the shared library has moved, or in the case of Boost's generated libraries, it compiled in a relative path rather than an absolute path. I ran into this issue before with CUDA's dynamic library (I was flabbergasted that it didn't automatically store an absolute path when I passed an absolute path to the linker), and my half-kludge of a fix was to add the environmental variable DYLD_LIBRARY_PATH=:/usr/local/cuda/lib.

Well, it happened again today; Boost was failing with dyld not loaded, image not found. I hit upon the google query that brought up an Apple developer page that explained the issue well.

I invented a general solution, since this will probably happen again to me. If you have shared libraries that don't use an absolute path (or don't use the right path), cd to the directory in which they reside and run the following command:

for f in *.dylib; do sudo install_name_tool -id $(pwd)/$f $f; done

Ta-da!

Ah, keyboard shortcuts.

I'm too used to using vim, where Ctrl-w deletes the preceding word in insert mode. (The same shortcut works in the terminal and other places, as well.) The downside is that when I'm on a Windows machine, i.e. using Firefox at the library, trying to delete a word reflexively closes the window in which I'm typing. This has happened three times already and almost happened once while typing in this post.

Abbreviating Russian names in BibTeX

Some of my research involves a Russian-developed method called Quasidiffusion, so naturally the appropriate works need to be properly cited. BibTeX, however, was mangling the correctly abbreviated Russian names.

"Дмитрий Ю. Анистратов" is transliterated as "Dmitriy Yu. Anistratov." Note that the middle initial is one character in the Cyrillic alphabet, but it becomes two Latin characters. Even when I enclosed the middle initial with a set of braces (which in most cases forces it to represent the text exactly), BibTeX truncated it to one Latin character (ANISTRATOV, D. Y.), which is incorrect.

The solution, gleaned from this document, is to use {\relax Yu}. as the middle name. Now, the author entry Anistratov, Dmitriy {\relax Yu}. is correctly abbreviated as "ANISTRATOV, D. YU."

Trilinos, MPI, Boost on my mac

Using CMake as the foundation of my code's build system has really made my life easier. When I add a compiler to /usr/local/bin, even though my PATH is set to prefer binaries in that directory, my old build directories will still compile with the original Apple-installed gcc in /usr/bin, all because CMake uses absolute paths for everything set at configuration time.

Experience has also shown me that, when installing software packages, it's best to put them in different subdirectories of /usr/local . CMake also makes using this structure much less of a hassle, since it's easy to specify the "-L" directory for each package.

Anyway, today I downloaded the pre-compiled GCC 4.4 package from the mac HPC site and installed it to its default location of /usr/local . Next, after setting CXX=/usr/local/bin/g++ and similarly with CC, FC, and F77, I downloaded and installed the latest version of Open MPI, with --prefix=/usr/local.

Finally, I installed Trilinos to a new path at /usr/local/trilinos-9.0.2-parallel/:

./configure
--cache-file=config.cache \
--prefix=$INSTALLDIR \
--enable-mpi \
--with-mpi-compilers \
CFLAGS="-O3 -ftree-vectorize" \
CXXFLAGS="-O3 -ftree-vectorize" \
FFLAGS="-O2 -ftree-vectorize" \
--with-libs="-framework vecLib" \
--with-ldflags="-Wl,-multiply_defined -Wl,suppress" \
--enable-epetra \
--enable-aztecoo \
--enable-amesos \
--enable-examples \
--enable-didasko \
--enable-teuchos \
--enable-triutils \
--enable-galeri

Now, after making sure that mpirun points to the correct version in /usr/local, I can link against the Trilinos libraries when using the auto-vectorizing GCC 4.4.

I also compiled Boost with the new GCC into its own directory (bjam file at ${HOME}/user-config.jam:

using mpi ; using gcc : 4.4 ; 
using python : 2.5 : /opt/local/bin :
 /opt/local/include/python2.5 : /opt/local/lib : ;

) using bjam --with-mpi --with-python --with-math --with-serialization --with-test --toolset=darwin install.

For some reason, though, even when I link against the libraries using an absolute path (/usr/local/bin/g++ -O3 -DNDEBUG -Wl,-search_paths_first -headerpad_max_install_names -fPIC CMakeFiles/basicTest.dir/basicTest.cpp.o -o ../bin/basicTest ../lib/libimclib.a ../lib/libmclib.a /usr/local/lib/libmpi_cxx.dylib /usr/local/lib/libmpi.dylib /usr/local/lib/libopen-rte.dylib /usr/local/lib/libopen-pal.dylib /usr/lib/libutil.dylib /usr/local/boost-gcc44/lib/libboost_mpi-xgcc44-mt-1_39.dylib /usr/local/boost-gcc44/lib/libboost_serialization-xgcc44-mt-1_39.dylib ), it doesn't store the absolute path to the shared library in the binary. That is to say, otool -L basicTest spits out:

 /usr/local/lib/libmpi_cxx.0.dylib (compatibility version 1.0.0, current version 1.0.0)
 /usr/local/lib/libmpi.0.dylib (compatibility version 1.0.0, current version 1.0.0)
 /usr/local/lib/libopen-rte.0.dylib (compatibility version 1.0.0, current version 1.0.0)
 /usr/local/lib/libopen-pal.0.dylib (compatibility version 1.0.0, current version 1.0.0)
 /usr/lib/libutil.dylib (compatibility version 1.0.0, current version 1.0.0)
 libboost_mpi-xgcc44-mt-1_39.dylib (compatibility version 0.0.0, current version 0.0.0)
 libboost_serialization-xgcc44-mt-1_39.dylib (compatibility version 0.0.0, current version 0.0.0)
 /usr/local/lib/libstdc++.6.dylib (compatibility version 7.0.0, current version 7.12.0)
 /usr/local/lib/libgcc_s.1.dylib (compatibility version 1.0.0, current version 1.0.0)
 /usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 111.1.4)

which means that the program fails when trying to load, because the library search paths don't include /usr/local/boost-gcc44/lib. The solution I came up with, after reading through documentation for Apple's ld and finding nothing to help, was to just symlink Boost's libraries into /usr/local/lib. Fortunately, it doesn't conflict with my gcc40 binaries because of Boost's name mangling scheme.

So now, by specifying absolute paths to the compiler, to mpirun, and to the libraries, I can compile with either Apple's compilers or with my own. Woot.

It's crazy how much of this build process hacking I've come to understand over the last couple of years. When I was at Oak Ridge, I don't think I even understood what a shared library was, or why my linker would spit out errors when I put the wrong path to a library.

Checking vectors of doubles with Google Test framework

The Google C++ testing framework makes writing unit tests (with helpful output) much easier. One of the capabilities it lacks is to verify that two containers of doubles (usually std::vector) are equal in each element and in their size.

It will work on any iterable container of doubles: lists, vectors, Blitz::TinyVectors, etc.
It will also check that the lengths are the same and verify that one is not shorter than the other, and can do it without assertions that cause the unit test to prematurely abort.

#define EXPECT_ITERABLE_DOUBLE_EQ( TYPE, ref, target) \
{ \
const TYPE& _ref(ref); \
const TYPE& _target(target); \
TYPE::const_iterator tarIter   = _target.begin(); \
TYPE::const_iterator refIter = _ref.begin(); \
unsigned int i = 0; \
while(refIter != _ref.end()) { \
    if ( tarIter == _target.end() ) { \
        ADD_FAILURE() << #target \
            " has a smaller length than " #ref ; \
        break; \
    } \
    EXPECT_DOUBLE_EQ(* refIter, * tarIter) \
        << "Vectors " #ref  " (refIter) " \
           "and " #target " (tarIter) " \
           "differ at index " << i; \
    ++refIter; ++tarIter; ++i; \
} \
EXPECT_TRUE( tarIter == _target.end() ) \
    << #ref " has a smaller length than " \
       #target ; \
}

Call it with something like:

 EXPECT_ITERABLE_DOUBLE_EQ(std::vector, expectedTemp, radTemperature);

and if there's a problem, you'll get an error like:

tRadiationController.cpp:163: Failure
Value of: * refIter
  Actual: 1.0999999999999992
Expected: * tarIter
Which is: 666
Vectors expectedTemp (refIter) and radTemperature (tarIter) differ at index 1

EDIT 7/7/09: Modified to allow passing the results of function calls into the arguments. This would unquestionably be better written as a templated function rather than a macro, except that we would not have the correct line number/file name, and our error message would not have the variable names.

CMake, Python, and MacPorts

Over the last few days, I've been learning the CMake build system. There's a lot that I like about it—it's certainly more scalable than SCons, which I've been using until now. In about a day of fooling around from scratch, I converted one of my main projects (MCGeometry) to CMake.

However, some stuff that should be trivial is tricky (the most frustrating thing is trying to copy files to a build directory -- it installs stuff to the prefix path without a problem, but intermediate copies don't work well at all). I've gone over the fooling around needed to get boost.python to work with MacPorts, but one would expect CMake (since its configuration allows the manual selection of most compilers, libraries, etc.) to handle it without a problem. Sadly, though, when using the FindPythonLibs module and FindSWIG and UseSWIG, it defaults to linking against the Python framework found in /System/Library. I spent probably a few hours googling for a solution, but finally today I found one that works without any hacks at all.

The key is that on Mac OS X systems, rather than looking for headers and libraries, FindPython prefers to search for the Python framework. There are several locations it checks in some preferential order. All that needs to be done is to symlink the Python framework installed by MacPorts to one of the higher-preference locations:

ln -s /opt/local/Library/Frameworks/Python.framework /Library/Frameworks/

Ta-da, now running otool -L on the SWIG modules that CMake generates properly shows that they link against the MacPorts framework, and the right include and library paths are used.

Boost mpi and python with macports

On my system, I installed MacPorts (formerly DarwinPorts) so that I could get matplotlib and numpy without any extra effort. In the process, it installed its own version of Python in /opt/local instead of the default mac one at /System/Library/Frameworks/Python.framework/Versions/2.5/. The former is now the only Python version I use on my computer. It hadn't posed a problem until I tried to install boostmpi (see previous post on MPI with C++ and Python).

Well, after I told it to use the correct I got the unfortunate error:
"Fatal Python error: Interpreter not initialized (version mismatch?)" A little searching showed that the problem was that the libraries (specifically, the libboost_python dynamic library) was linking against the dynamic library found in /System/Library rather than the one set up in my environment and everywhere else. Using otool -L on the dynamic lib confirmed this, and I used the same tool while trying to find a solution (rather than using the Python failure to check).

On the Mac, installing Boost properly (with any of the libraries that need compiling) requires passing the --toolset=darwin option to bjam. Grepping the Boost source for "/System/Library" revealed that location to be hardwired as the default Python location if the darwin toolset is used. The solution is to specify in the user-config.jam file exactly what version of Python to use and where. My final install script (echo, bjam, and rm are the three commands):
echo "using mpi ; using python : 2.5 : /opt/local/bin : /opt/local/include/python2.5 : /opt/local/lib : ;" > ${HOME}/user-config.jam
bjam --with-mpi --with-python --with-math --with-serialization --toolset=darwin install --prefix=$INSTALLDIR
rm ${HOME}/user-config.jam

Using otool -L on the new libboost_python-xgcc40-mt-1_39.dylib verified that it linked against the correct Python, and boostmpi ran correctly after installing it.

Python MPI, with C++

The front end to my experimental code is in Python (with C++ running all the important stuff, and SWIG coupling them), which is how I want it to stay. Sooner or later my problems are going to outgrow my computer, and I'll need to run them on a cluster in parallel using MPI. It wasn't clear whether this was possible, or how it was to be done. But Python at its heart is C, and the shared objects it runs with are just code. I knew that, for example, Trilinos was somehow able to use both MPI and Python; so it's certainly possible.

The first thing I did was to search for implementations where people just used MPI with Python (without worrying about the C++ side). There are several such implementations, including the apparently stillborn pyMPI, the numpy-dependent Pypar, the basic mpi4py, and Python bindings to MPI through Boost::Python and Boost.MPI. Boost.MPI looked like the best implementation of them all, but unfortunately I could not for the life of me figure out where it was in the distribution or how to install it. (It turns out the version I was using, 1.36, didn't have it, but even in the latest 1.39 I can't see how to install the Python module.)

So the first task I did was to install mpi4py, which was trivial, and to write a quick program. mpi4py is a laughably thin front-end, so the code is not at all Pythonic.
from mpi4py import MPI
comm = MPI.COMM_WORLD
print "Hi, I am proc %d of %d" % (mpi.rank + 1, mpi.size)
Then, to run it, mpirun -np 2 python run.py

That worked fine. Surprising to me was how easily the MPI calls in C++ worked with this. I created a small class, Thing, with a member function Thing::sendRecv() that did some MPI_Send and MPI_Recv calls (depending on whether the processor was rank 0 or not), and return. I put a quick SWIG wrapper around it, compiled it to generate _thing.so and thing.py, and changed my Python code to:
from mpi4py import MPI
import thing
comm = MPI.COMM_WORLD
print "Hi, I am proc %d of %d" % (comm.Get_rank(), comm.Get_size())
t = thing.Thing()
t.sendRecv()


Again, mpirun -np 2 python run.py, and it worked! The subroutine did its thing and returned with no problem. It looks like using both MPI and Python with scientific computing is good to go.

After this was working, I decided to try to use the "standalone" boost.mpi, but that was opening up a whole other can of worms. My excursion into the land of otool and bjam will be detailed soon. Suffice it to say that, once I got boost.mpi compiled and linked against the right version of Python, the same code (except the MPI in Python looked much better).

I surmise that the use of any Python MPI library can be rendered unnecessary as long as your C++ code calls MPI_Init and MPI_Finalize in the right places.

MCNP criticality search Python script

One of the occasional tasks of a nuclear engineering student is to perform a criticality search using some parameter. That is, by varying (for example) the radius r of a uranium sphere, at some point the neutron production (from fission) and loss (from absorption and leakage) terms exactly balance, and the reactor is critical. The issue is complicated a little when using a Monte Carlo simulation, since the k-eigenvalue (unity when critical) is only reported with a certain probability or tolerance; running more particles will decrease the uncertainty.

Undergraduates, being unthinking animals, usually try this arduous task by hand, laboriously changing the parameter in an input file, running it, and repeating. As a graduate student, such tediousness is beneath me, as I would rather spend an equal amount of time writing a script to do the task for me. (Experience dictates that most of the time, some error or revision will require re-doing the task at least once, which means that a script will often save time by a factor of two or more. Besides, scripts are reusable and an educational experience.)

These Python scripts (gzipped archive) will do just that. It's a fairly simple binary search (no Newton–Krylov iterations) partly because of the complicating factor of the uncertainty in each answer. When criticality is reported to within the experimental standard deviation, the script increases the number of particles until such a parameter is determined that results in criticality within the user's specified tolerance. The script will run MCNP, process the output to determine k-effective, and iterate to your heart's content. If you have matplotlib installed, it will also produce a plot of how k varies with the parameter:

The MCNP search script is by no means perfect, and it assumes you have mcnp5 in your path. You have to know some Python to create a subclass of McnpSolver for your particular problem.

Included are scripts for finding the critical parameters of HEU cubes, cylinders, hemispheres, etc. There's also a script for generating a spherical shell volume tally to calculate the radial flux of the converged critical sphere.

Future updates of the MCNP criticality script will be posted on my projects page, although you should not expect there to be any. If you use this script, please cite its use if you report or present the results.

En dash

Today, having noticed the behavior in my literature review, I learned that en dashes are the appropriate punctuation to use between the names of joint authors. For example, "Newton-Krylov methods" should be rewritten as "Newton–Krylov methods."

See this google books result.

Overwriting lines in terminal output

For weeks I've been wondering, off and on, if there existed an easy way to write information to the terminal output and be able to change it later, to overwrite or overtype or erase or something. I'd seen how curl will update the fraction downloaded and download rate at the end of the line rather than print a new line for every status change; such a feature would be very nice my code, where I could just use one line of output to write information about the time step. I'd seen the "curses" library, but that looked ridiculously complex for what I was asking. (It creates a whole GUI in the terminal rather than allowing simple operations.)

Well, after finally hitting the right keywords in a google search, I came across Python with ANSI escape codes (more ANSI escape sequences here). This is just what I want. All you have to do is to send the '\r' carriage return sequence to move to the beginning of the line, and then another ANSI sequence can erase the line. That behavior is actually just what I would expect based on the descriptions of the escape sequences: '\n', the "new line" character, should create a new line; '\r', which is called the 'carriage return', actually *does* in this contact emulate the behavior of a typewriter carriage returning to the rightmost position. Mostly due to the cross-platform confusion of the "end-of-line" character business, I'd never really thought this was the way it would work. The long and short of it, in a few lines of python code:

from time import sleep
from sys import stdout
for i in range(10):
    stdout.write('\r\x1b[K')
    stdout.write('Time step %03d/%03d  ' % (i+1, 10))
    stdout.flush()

The '\r\x1b[K' moves the cursor to the beginning of the line the following ANSI escape sequence clears the line completely. The next line actually writes the code. The trailing spaces are in case some sort of output gets print during the time step: it will go at the end of the line, and assuming it has a newline of its own, the time step at which the output occurs will be preserved on the screen, and the next line will have the actual updating time step information.

Lessons learned of the evening

To make a long story short, or more precisely, to summarize an hour of frustration: when using SWIG with Python, watch that you actually instantiate an object by calling its constructor, because apparently you can change the values of the class itself without creating an object from it. Then, when you later want to pass an instantiated object to another C++ method, it won't fail with an inscrutable error. You probably don't want to read the rest of this post, but it's preserved for posterity.

I have a struct of data, SourceDataUniform, that I pass to a class called RadiationController. My original code was:

newSource = ImcTransport.SourceDataUniform
newSource.emissionRate = 2.0
newSource.cellIndex = 1
self.__radController.addUniformSource(src)

but this failed with TypeError: in method 'RadiationControllerT_addUniformSource', argument 2 of type 'iMc::SourceDataUniform'. Yet I could print the data in the object, and it said it was of type <class 'ImcTransport.SourceDataUniform'>.

I knew I'd gotten the exact same thing to work with a different class, so I wasn't sure if the problem was that I was using a struct, or what.

I spent a good while looking through the SWIG documentation on proxy classes, structs, the ".this" method. Eventually, noting that the similar case that actually worked gave a repr of <ImcPhysics.OpacityModelConstant; proxy of <Swig Object of type 'iMc::OpacityModelConstant *' at 0x305670> >, and realizing that the "this" object was only set in the constructor for the SourceDataUniform Python wrapper, I discovered that I wasn't instantiating an object. Adding parentheses

newSource = ImcTransport.SourceDataUniform()

solved the problem.

Unit tests again

Life is funny sometimes. It's occasionally the case that one thinks some almost-trivial task is totally unnecessary, yet the "best practice" is to do it anyway&emdash;and the task turns out to catch or prevent a subtle but insidious error.

This evening I started writing unit tests for my code, and the first one on the list was the (almost) dead-simple "check direction vector" routine, which is overloaded on Blitz tinyvector length. The one-D case has to ensure that the cosine is bounded by one. Simple enough, right?

inline bool checkDirectionVector(
        const blitz::TinyVector<double, 1>& omega)
{
    return (abs(omega[0]) <= 1.0);
}

Well, in addition to testing that it called a value of 1.0 and -0.5 allowed, I tested whether -1.1 was disallowed. I was flabbergasted when I got the message

  /---- Unit test tBlitzStuff failed
  |   in <build/imclib/tests/tBlitzStuff.cpp> on line 32
  |   !checkDirectionVector(mu)
  +-------------------------------

which corresponded to the code

    mu = -1.1;
    TESTER_CHECKFORPASS(!checkDirectionVector(mu));

What the heck? Then I realized the mistake. My compiler (gcc 4.0) was silently downcasting a double to an integer, truncating -1.1 to -1, and then upcasting the result to the double -1.0, which satisfies the condition.

The solution is simple, to just use std::fabs(omega[0]) instead. But seriously? I never would have expected a failure on that unit test of all places. These things aren't as worthless as one might think.

Zoidberg, on the economy

Once again, the conservative, sandwich-heavy portfolio pays off for the hungry investor!

Ok, I like Python now

Whoever said that Perl had to be the language in which to write super-complicated commands all on one line? This Python command will turn a list of lists into a string amenable to writing to a file as MATLAB input:
result = '[' + "\n ".join([" ".join(["%6g" % val for val in line]) for line in array]) + ']'

EDIT: using generator comprehensions rather than list comprehensions for greater efficiency (thanks, Mike):result = '[' + "\n ".join(" ".join("%6g" % val for val in line) for line in array) + ']'

More frustration (Vim) resolved

Ugh, I spent the entire morning trying to fix a problem with Vim. I remember at one point getting frustrated that it refused to handle indentations of comments (which start with the hash or pound sign, #) correctly. The solution I got was to use someone else's indent file. Well, that indent file doesn't handle formatting lists extended over more than one line. But when I went back to the original, which worked for that, it couldn't handle the comments. I spent probably three hours searching for a solution, looking through $VIM/runtime/ trying to find out what the differences were between the different indent files. Furthermore, although the 'sh' file type displayed this behavior, the 'perl' file type did not!

Well, the long story was made short when I realized that the reason it forced all lines with a '#' to the beginning: it thought they were C preprocessor directives, which belong in column 1. I had unknowingly set an option, smartindent, which caused this behavior. Somehow, and I still haven't figured this out, the options in the Perl file mitigated that behavior.

Moral of this story: avoid setting options in vim that you aren't REALLY sure you understand. And with vim, this probably means you won't be setting many.

Template specialization for std::swap

One of the algorithms in C++'s STL is swap, a templated function which takes two arguments of the same class and switches them. By default, it creates a new object and performs two copies, but for STL containers (such as std::list, which I was using) it has its own template specialization to only copy their internal pointers and not any actual data. This makes, by default, std::swap way more efficient on an STL container than everything else.

I have a class Census inside a namespace iMc, and it has-a std::list to store a bunch of particles. I'd like to call std::swap(iMc::Census& a, iMc::Census& b) but not copy all the particles inside the list. The solution was obvious and allowed: create a new specialization for my own class that calls  std::swap on the std::list :

namespace std {
template<>
inline void swap<iMc::Census>(iMc::Census& a, iMc::Census& b)
{
std::swap(a._particles, b._particles);
}
} // end namespace std

Well, because the _particles list was private, I had to make this function a friend of iMc::Census. Easy, right? Wrong. It turns out everything on the internet; that talks about templated friend functions is a lie! No matter what I did, I received a bunch of weird errors.

namespace iMc {
class Census {
private:
template<typename T> friend void std::swap(T& a, T& b);
...
};
}

gives

build/imclib/Census.hpp:64: error: 'void std::swap(iMc::Census&, iMc::Census&)' should have been declared inside 'std'

and

friend template<> void std::swap<iMc::Census>(iMc::Census& a, iMc::Census& b);

gives

build/imclib/Census.hpp:64: error: expected unqualified-id before 'template'

Well, it turns out that the correct answer is to have

namespace iMc {
class Census {
private:
template<typename T> friend void std::swap(T& a, T& b);
...
};
}

Who would have guessed? I even checked Stroustroup's reference book, and it didn't look like it had syntax for templated friend functions in there. Hopefully no one else has to fool around and waste an hour of their time trying to do this.

CC license

I'm LaTeXing my class notes, and I plan on distributing them in some sense. Here's the Creative Commons license in LaTeX that I came up with (which puts CC in a circle and uses the hyperref package to hyperlink it):

\newcommand{\cclicense}{%
  \smash{\raisebox{-0.45ex}{%
  \setlength{\unitlength}{1em}%
  \begin{picture}(1,1)%
    \put(0.5,0.5){\circle{1}}
    \put(0.5,0.5){\hbox to 0pt{%
      \hss\raisebox{-.45ex}{\tiny\textsf{CC}}\hss}}%
  \end{picture}%
  }}%
  \hskip -1em%
  \href{http://creativecommons.org/licenses/by-sa/3.0/}%
  {\ \hskip 1em \textsf{BY-SA}}%
}

Update: thanks are due to Mike for pointing out my use of an externally defined term, which I thought was in base/latex.ltx but was actually defined elsewhere.

MCNP syntax highlighting in vim

This syntax highlighter for MCNP in VIM is due to Giacomo Grasso. The file on his page is inaccessible other than through the white paper, so I have taken the liberty of extracting it and putting into a file available for download. Put it inside the ~/.vim/syntax folder.

There's no easy sure-fire way of detecting an MCNP input file, but if you create a file at ~/.vim/scripts.vim with the contents

if did_filetype() " filetype already set..
  finish  " ..don't do these checks
endif
if getline(2) =~ '^C.*\<mcnp\>\c'
  setfiletype imcnp
endif

it will load the MCNP syntax coloring if the second line of your input deck is a comment (the first letter of the line is a 'c') with the word MCNP in it somewhere.

Return of the "pickup notice"

I have gotten one or two of these in the past. This time it was a lot less mysterious. Even though it had a legitimate non-bulk stamp on it, it was labeled with "enroll at www.feinfo.org" and had a sticker that told me to "CHOOSE [MY] GIFT." Among my choices were a "condo stay," an "iPod shuffle," and a "visa gift card." This one was also labeled with a copyright by "NETWORK DIRECT, INC" which is one of these bogus advertising companies. I've heard about these "just show up to a meeting and you'll get a prize" before. They've been around for decades, and you have to resist a salesman's insistent pushes for an hour or more. That doesn't sound pleasant, and I imagine that if it were too good to be true, it probably would be.

Vim/TeXShop integration

I like to edit my LaTeX documents in vim (with the VIM-LaTeX package), but for various reasons, I like to use TeXShop to actually compile and view the document. In TeXShop, make sure that the "Configure for External Editor" button is checked, and then quit it. Once you add the following lines to your ~/.vim/ftplugin/tex.vim file, you just have to hit Cmd-R for the document to automatically open in TeXShop, and compile. You can use everything else just as you would expect. It really works great.

" Run LaTeX through TexShop
function! SRJ_runLatex()
   if &ft != 'tex'
      echo "calling srj_runLatex from a non-tex file"
      return ''
   end

   "write the file
   :w

   let thePath = getcwd() . '/'. expand("%")

   let execString = 'osascript -e "tell app \"TeXShop\"" -e "set theDoc to open ((POSIX file \"'.thePath.'\") as alias)" -e "tell theDoc to latexinteractive" -e "end tell"'
   exec 'silent! !'.execString
   return ''
endfunction
no  <expr> <D-r> SRJ_runLatex()
vn  <expr> <D-r> SRJ_runLatex()
ino <expr> <D-r> SRJ_runLatex()

The <D-r> means that hitting Command-R will save and compile it. You can change it in the syntax of any of the other vim command shortcuts.

EDIT 3/24/2009: To make it compatible with 10.4, since apparently the "open" command does not return the opened document as it does in 10.5, change the osascript line to:

let execString = 'osascript -e "tell app \"TeXShop\"" -e "set theDoc to open ((POSIX file \"'.thePath.'\") as alias)" -e "try" -e "tell theDoc to latexinteractive" -e "on error" -e "set theDoc to front document" -e "tell theDoc to latexinteractive" -e "end try" -e "end tell"'

UPDATE 9/25/2011: I've switched from using TeXShop to Skim. See my new Vim-LaTeX/Skim scripts here

Insolation pt. 2

Incidentally, the plot I posted yesterday was only of the intensity averaged over the daylight hours. Multiplying by the hours of daylight (generated with the same program) should give the total amount of solar energy deposition per unit area in a day.

UPDATE This actually seems to give the wrong values, according to here and here. Maybe the NASA calculations are weighted differently than they say, or I'm misunderstanding what they mean. Either way, ignore this plot and use the values shown by those links.

Insolation

Today's topic is not "insulation," as in "my apartment doesn't have enough insulation to keep me warm in this frigid weather." Rather, insolation is how much sunlight hits the earth's surface.
Ever since moving up north and noting how low the sun is in the sky and how short are the days during winter, I've wondered about how much actual energy gets dumped onto a unit area on the surface of the earth.

After googling all sorts of permutations on "day integrated solar flux" and similar terms, I stumbled upon the more accurate term of insolation (which, after the fact, I remembered having heard through science bowl). Searching around some more resulted in this factual page, and further searching led me to a NASA page with all sorts of useful information. Based on values from one of NASA's programs, which takes into account most of the important stuff, here is a plot of the solar intensity as a function of the days of the year:

It's no wonder it feels dark and cold here in the winter. A similar plot of the average solar angle shows how the sun is always lower in the sky here, too.

Of course, the plot does not include the extra atmospheric attenuation from the sun coming in at lower angles. (The same thickness of clouds will require sunlight to travel further as a result.) That gives Ann Arbor and the other northern places an even more severe disadvantage.

Incidentally, the MATLAB code to display the months (where the x axis data is stored as days) is:
monthlengths = [0 31 28 31 30 31 30 31 31 30 31 30];
months = {'J' 'F' 'M' 'A' 'M' 'J' 'J' 'A' 'S' 'O' 'N' 'D'};
set(gca,'XTick',cumsum(monthlengths));
set(gca,'XTickLabel',months);
and I created the graphic with the same script I posted about earlier.