Using Boost Spirit 2.1+ to evaluate constant arithmetic expressions

I used Boost Program_options for my research code's input parsing requirements. As expected, it works well. For floating point quantities, I've found having to specify 4.18879020478639 for 4π/3 is a bit cumbersome. When specifying arguments on the command line, I could use any one of a number of hacks to perform floating point arithmetic in bash. But that didn't help me when taking input from response files. I wanted an arithmetic expression utility that could process a C-like grammar. Long double accuracy was a requirement. I looked around a bit. GNU libmatheval and muParser appeared to be good fits but looked to be double precision-only. Though straightforward, I decided against combining the Shunting-yard algorithm with some in-place RPN logic because I wasn't enthusiastic about either finding or writing an appropriate (and admittedly toy simple) lexer.

Already having a Boost pre-requisite in my code, I decided to try out Boost Spirit since it could provide a header-only, precision-agnostic solution. Spirit has gone through a lot of significant revisions and many articles seemed to be written for the older Spirit Classic. The latest official reference manual gave several excellent examples that almost fit. I adapted Spirit's calc3.cpp and roman.cpp sample grammars to accomplish what I wanted. I added ** as an operator for exponentiation and pi as a built-in constant. I also added a smattering of special functions (e.g. exp). The task required about 150 lines of Spirit-based code. Further extensions should be fairly straightforward.
My sample code does the following:

1. Declares X-macros for the unary and binary functions built into the grammar.
2. Declares Boost Phoenix-ready functors for the unary and binary functions.
3. Declares expression::grammar consisting of the rules expression, term, factor, and primary.
4. Declares a helper called parse to hide needing boost::spirit::ascii::space when performing whitespace-agnostic parsing.
5. Declares a small test macro called EXPRTEST to reduce test case boilerplate.
6. Runs a collection of tests to ensure the parser is behaving as expected.

Compared to the samples on which I built, the biggest changes to the grammar were adding primary and reworking factor. I could have shortened up some details with a few using boost::spirit::qi; declarations, but I wanted to keep track of exactly which parts of Spirit I used. Feedback and suggestions are most welcome.

Updated 9 Jan 2012: During a second pass motivated by reducing compilation times, I removed the X-macro nastiness and replaced it with some qi::symbols-based logic. The code in the anonymous namespace could likely be removed if I could figure out the right boost::phoenix::bind hoodoo to use to replace lazy_ufunc_ and lazy_bfunc_ (any suggestions?). I have not yet measured the difference in either memory footprint or runtime performance. Compilation time and memory is greatly reduced for both GNU and Intel compilers.

Updated 31 Mar 2013: I have moved the code into a gist to facilitate both maintenance and gathering feedback. I have added the constant e and slightly reordered the primary grammar per the comments.

Updated 17-19 Sept 2013: I have added an MPL2 header to the gist. In the comments is a suggestion for how to use standalone MPL2 logic building atop this grammar from my project Suzerain to incorporate such parsing into your own, possibly non-MPL2, source tree.

Adding enthought.traits.api.HasTraits as an ancestor using Python metaclasses

Today, while banging my head against a wall trying to get a SWIG-generated Python module to play nice with Python Traits, I cooked a way to use Python metaclasses to modify a class to add enthought.traits.has_traits.HasTraits as a base class:

from enthought.traits.api import HasTraits
from enthought.traits.has_traits import MetaHasTraits

# See http://onlamp.com/lpt/a/3388 for an intro to Python metaclasses and
# http://code.activestate.com/recipes/204197-solving-the-metaclass-conflict/
# for what can go wrong here
class AddHasTraitsAncestor(MetaHasTraits):
    def __new__(cls, name, bases, dct):
        return MetaHasTraits.__new__(cls, name, cls.modify(bases), dct)

    def __init__(cls, name, bases, dct):
        super(AddHasTraitsAncestor, cls).__init__(name, cls.modify(bases), dct)

    @staticmethod
    def modify(bases):
        if (bases == (object,)):
            return (HasTraits,)
        else:
            new_bases = list(bases)
            new_bases.insert(0, HasTraits)
            return tuple(new_bases)

Adding the line __metaclass__ = AddHasTraitsAncestor to a class definition is enough, in some cases, for the class to now have HasTraits as a base class. Most of the frustrating bits were before I realized HasTraits has its own metaclass MetaHasTraits which AddHasTraitsAncestor had to subclass before the approach worked.

"Uhhh...?" you say. That puzzled look asking "Why didn't you just add the base class yourself?". I used this approach because SWIG was generating the Python class definitions and I was abusing SWIG's pythoncode feature to tweak the SWIG wrappers to be Traits-ready a la

%extend someCPlusPlusClass {

    %pythoncode %{
        __metaclass__ = AddHasTraitsAncestor
    %}

}

Turns out it works and my SWIG wrappers are now Traits-friendly. I'm still working on getting Traits to see the SWIG-based members, however. Suggestions much appreciated.

Using gnuplot from the command line

gnuplot is excellent for both quick 'n' dirty and lengthy plotting projects. Still, for truly quick work of just plottin' some column against some other column, I've always felt like gnuplot demands too many keystrokes:

[1000 rhys@tiresias 16]$ gnuplot

 G N U P L O T
 Version 4.2 patchlevel 6
 last modified Sep 2009
 System: Linux 2.6.32-30-generic

 Copyright (C) 1986 - 1993, 1998, 2004, 2007 - 2009
 Thomas Williams, Colin Kelley and many others

 Type `help` to access the on-line reference manual.
 The gnuplot FAQ is available from http://www.gnuplot.info/faq/

 Send bug reports and suggestions to 


Terminal type set to 'x11'
gnuplot> plot 'restart0.mean' using 2:17
gnuplot> exit

Of course, one can turn this process into a one-liner at the cost of quoting:

gnuplot -e "plot 'restart0.mean' using 2:17"

Call me picky, but the signal-to-noise ratio in that one-liner still feels too high.

After some quoting experimentation using printf's %q and %b format specifiers, the following little function definition

gplot () {
   local fileexpr=$(printf "'%q'" $1)
   shift
   local plotexpr=$(printf '%b' "plot ${fileexpr} $*")
   gnuplot -persist -raise -e "$plotexpr"
}

makes

gplot restart0.mean using 2:17

behave like the previous two examples. This function happily hides the file name quoting requirements and remains tab-completion-friendly.

Providing computed column values like gnuplot's using 1:$(cos($2)) on the command line still requires thought, but the straightforward tasks work thoughtlessly. Suggestions for improvements definitely welcome in the comments. One glaring bug is that passing a pipe expression as the file bombs.

Update 20110504: I've decided that I much prefer creating a gplot script containing

#!/bin/bash
file=$1
shift
gnuplot -persist -raise <(cat <<-GNUPLOTSCRIPT
plot '$file' $*
GNUPLOTSCRIPT
)

as it is much simpler to work through the quoting and far more extensible.

Update 20111031: Globbing, labelling, and a handful of other enhancements:

#!/bin/bash

# Fail on first error
set -e

# Create temporary files to hold gnuplot script
tmp1=`mktemp`
tmp2=`mktemp`
trap "rm -f $tmp1 $tmp2" EXIT

# Process (and then remove) command line options
# Build up any post-terminal processing in tmp2 for later use
autotitle=
title="`date +%r`: gplot $*"
terminal="set terminal x11 title '$title' enhanced persist"
cmd=plot
forexpr=
raise=
showhelp=
while getopts "3cf:ghils:re:t:x:y:z:F:SX:Y:Z:" opt
do
  case $opt in
    3) cmd=splot
       ;;
    c) autotitle=true
       ;;
    e) terminal='set term postscript eps enhanced color dashed rounded "Arial,14"'
       pause='' # Disable interactive on EPS output
       echo set output \"$OPTARG\" >> $tmp2
       ;;
    f) forexpr=" for [$OPTARG] "
       ;;
    g) echo "set grid" >> $tmp2
       ;;
    h) showhelp=0
       ;;
    i) pause='pause -1 "Plot interactive until Enter pressed.\nHit h in plot window for help.\n"'
       ;;
    l) echo "set logscale y" >> $tmp2
       ;;
    r) raise="-raise"
       ;;
    s) savescript="$OPTARG"
       ;;
    t) echo set title \"$OPTARG\" >> $tmp2
       ;;
    x) echo set xlabel \"$OPTARG\" >> $tmp2
       ;;
    y) echo set ylabel \"$OPTARG\" >> $tmp2
       ;;
    z) echo set zlabel \"$OPTARG\" >> $tmp2
       ;;
    F) pause="pause $OPTARG; replot; reread"
       ;;
    S) stdin=true
       ;;
    X) echo set xrange \[$OPTARG\] >> $tmp2
       ;;
    Y) echo set yrange \[$OPTARG\] >> $tmp2
       ;;
    Z) echo set zrange \[$OPTARG\] >> $tmp2
       ;;
  esac
done
shift $((OPTIND-1))

if [ x$showhelp != x ]; then
    cat <<-HERE
Usage: gplot [OPTION]... (FILE|EXTGLOB) GNUPLOTCMD...
Use gnuplot to plot one or more files directly from the command line.

  -3             Perform 3D plotting using gnuplot's splot command.
  -c             Populate the key using autotitling.
  -e FILE        Save an Encapsulated Postscript (eps) called FILE.
  -f FOREXPR     Prepend a 'for [FOREXPR]' to the plotting command.
  -g             Show grid lines.
  -h             Show this help message.
  -i             Interactive plotting mode.  Hit 'h' on plot for help.
  -l             Use logarithmic scale for y axis.
  -s FILE        Save the generated gnuplot as a script called FILE.
  -r             Have the window manager raise the plot window.
  -t TITLE       Set TITLE as the plot's title.
  -x XLABEL      Specify XLABEL as the x axis label.
  -y YLABEL      Specify YLABEL as the y axis label.
  -z ZLABEL      Specify ZLABEL as the z axis label.
  -F FREQUENCY   Replot the inputs every FREQUENCY seconds.
  -S             Prior to plotting, read auxililary gunplot from stdin.
  -X XLOW:XHIGH  Specify an explicit x axis range instead of autoscaling.
  -Y YLOW:YHIGH  Specify an explicit y axis range instead of autoscaling.
  -Z ZLOW:ZHIGH  Specify an explicit z axis range instead of autoscaling.

Examples (see gnuplot documentation for complete GNUPLOTCMD details):

  gplot -i foo.dat using 1:2 with linespoints
  gplot -s foo.gp -X 0:1 -Y 0:2 foo.dat using 1:2 with linespoints
  gplot -e foo.eps foo.dat using 1:2 with linespoints
  gplot -3 restart\*.dat using 1:2:3

On error, the failing gnuplot script is shown.
HERE
    exit $showhelp
fi

# Set terminal
echo "$terminal" >> $tmp1

# Slurp any settings built up during getops processing
cat $tmp2 >> $tmp1

# Obtain file(s) to plot from first argument using extended globbing
# Deliberately allow globbing to occur in this assignment
shopt -s extglob
declare -a files=($1)
shift

# Tweak autotitle based on options and incoming argument details
if [ "$autotitle" ]; then
    # Use columnheader option iff only one file provided
    if [ ${#files[@]} -eq 1 -a $cmd != splot ]; then
        echo 'set key autotitle columnheader' >> $tmp1
    else
        echo 'set key autotitle' >> $tmp1
    fi
else
    echo 'set key noautotitle' >> $tmp1
fi

# Possibly slurp standard input for further options
# FIXME Not working for 'echo foo | gplot'
if [ "$stdin" ]; then
    cat <&1 >> $tmp1
fi

# Build gnuplot script to plot possibly multiple files
declare -i count=0
for file in "${files[@]}"
do
    count+=1
    if [ $count -eq 1 ]; then
        echo "$cmd" "$forexpr" \'$file\' $* \\ >> $tmp1
    else
        echo "   "  , "$forexpr" \'$file\' $* \\ >> $tmp1
    fi
done
echo '' >> $tmp1

# If requested, save the plotting commands as an executable script
if [ "$savescript" ]; then
    echo '#!/usr/bin/env gnuplot' > "$savescript"
    cat "$tmp1" >> "$savescript"
    chmod a+x "$savescript"
fi

# If requested, add a pause command to the end of the script. Notice this was
# not added to the $savescript as it would be annoying in that context.
if [ "$pause" ]; then
    echo "$pause" >> $tmp1
fi

# Generate the plot
# Display the script when an error occurs to aid debugging
gnuplot $raise $tmp1 || cat $tmp1

Noise.

Following up on some quasi-1D success I'm ready to calculate some 3D fields. I seed them by taking my 1D solution and add noise. This is purportedly a terrible way to kick off turbulence as the fluctuations have to be massive to ultimately take and transition into a correct field. I've been warned...

Needing a way to monitor the fluctuation strengths and being in coefficient-space for a mixed Fourier/B-Spline discretization, I'm trying out the L^2 norm of the fluctuating signal. By Parseval's theorem along with some minor linear algebra, this turns into a sum over the complex-valued coefficient magnitudes:

Now, for the first case of the L^2 norm over all three directions I need the matrix M which computes the L^2 inner product of two functions given their B-spline coefficient expansions. The code to do this atop some old contributions to the GSL isn't bad, but testing that my banded derivative operator coefficients are correct is awful.

Happily, Mathematica 7 includes BSplineBasis which I can use to cook analytic results to test my implementation. Unhappily, BSplineBasis has issues computing derivatives and an awful amount of mucking around is required to get my analytical test cases:

In the end, this is all The Right Thing ™ to do but man it eats time. Recognizing this particular rabbit hole and how often such things ensnare my little mind, I sent the following email to my wonderful wife today:

In case you ever wonder what I do with my days....

1. I stare at things like this for far too long to work out the math.
2. Then I sit and stare at the math for far too long to figure out the computer implementation.
3. Then I stare at the computer implementation to test it and clean it up.
4. Then I use the implementation at a single line within my thesis research code.
5. Then I go talk to my advisor on Mondays and report that a single line now works this week that didn't work last week.

— Rhys