-fplan9-extensions

While reading 21st Century C, which shows some nice GCC extensions for stuffing anonymous struct members inside structs, I happened across -fplan9-extensions:

One can smell the homespun vtable lurking in the OOP use case...

Dos

About a week back, my wife gave birth to our second son Oxford C Thomas Ulerich.

First, pensive:

Second, at maximal yawn:

Driving a GSL Root Finder

The GNU Scientific Library one-dimensional root finding routines are handy, but I never remember how to drive them properly. Here is what I essentially always end up re-distilling from the examples:

int fsolver_solve(
        gsl_root_fsolver * const s,
        gsl_function     * const f,
        const size_t             maxiter,
        const double             epsabs,
        const double             epsrel,
        double           * const lower,
        double           * const upper,
        double           * const root)
{
    // Initialize fsolver on [lower, upper]
    gsl_error_handler_t * h = gsl_set_error_handler_off();    // Push handler
    int status = gsl_root_fsolver_set(s, f, *lower, *upper);
    status = (GSL_SUCCESS == status) ? GSL_CONTINUE : status;

    // Proceed until success, failure, or maximum iterations reached
    // On success, overwrite *location as described in API documentation.
    *root = GSL_NAN;
    for (size_t iter = 1; iter <= maxiter && status == GSL_CONTINUE; ++iter) {
        if (GSL_SUCCESS != (status = gsl_root_fsolver_iterate(s))) break;
        *lower = gsl_root_fsolver_x_lower(s);
        *upper = gsl_root_fsolver_x_upper(s);
        status = gsl_root_test_interval(*lower, *upper, epsabs, epsrel);
    }
    if (status == GSL_SUCCESS) {
        *root = gsl_root_fsolver_root(s);
    }
    gsl_set_error_handler(h);                                 // Pop handler
    return status;
}

Generating the Blasius boundary layer profile

I've been working on some boundary layer edge detection routines for Suzerain and wanted to use the Blasius boundary layer as a test case. The Blasius boundary layer produces a slick ODE that can only be solved numerically. Ganapol has a nice article on it and the more general Falkner—Skan equation. He also presents some nice tabulated data which I'd fit with an Akima, but I've found that (surprise surprise) the Akima natural end conditions disturb the second derivative of my fits sufficiently that it makes my test cases lousier than I'd like. Since Ganapol was the nicest tabulated Blasius data that I'd seen (meaning lots of digits, large similarity coordinate eta) and it didn't meet my needs, I set off to generate more (comparable number of digits, even larger similarity coordinate) using the GNU Scientific Library's ODEIV2 capabilities:

This code produces output matching Ganapol's everywhere to a minimum of eight digits. So though I'm reporting more than eight here, do not trust more than eight. The results are:

                   eta                       f                      fp                     fpp
0.0000000000000000e+00  0.0000000000000000e+00  0.0000000000000000e+00  3.3205733621519629e-01
2.0000000000000001e-01  6.6409997145970585e-03  6.6407792096250848e-02  3.3198383711462798e-01
4.0000000000000002e-01  2.6559884017994733e-02  1.3276416076102221e-01  3.3146984420145503e-01
6.0000000000000009e-01  5.9734637498038382e-02  1.9893725242221899e-01  3.3007912757429625e-01
8.0000000000000004e-01  1.0610822083904069e-01  2.6470913872311741e-01  3.2738927014925229e-01
1.0000000000000000e+00  1.6557172578927976e-01  3.2978003124966698e-01  3.2300711668694282e-01
1.2000000000000002e+00  2.3794871728892703e-01  3.9377610443395650e-01  3.1658919106111544e-01
1.4000000000000001e+00  3.2298157382953580e-01  4.5626176470513380e-01  3.0786539179016786e-01
1.6000000000000001e+00  4.2032076550162573e-01  5.1675678442261486e-01  2.9666346145571931e-01
1.8000000000000000e+00  5.2951803774381023e-01  5.7475814388945690e-01  2.8293101725975589e-01
2.0000000000000000e+00  6.5002436993528900e-01  6.2976573650238588e-01  2.6675154569727827e-01
2.2000000000000002e+00  7.8119333701057347e-01  6.8131037723602750e-01  2.4835091319037131e-01
2.4000000000000004e+00  9.2229012563454160e-01  7.2898193506257503e-01  2.2809176068668360e-01
2.6000000000000001e+00  1.0725059767878351e+00  7.7245502114856535e-01  2.0645462679942550e-01
2.8000000000000003e+00  1.2309773023143091e+00  8.1150962319910902e-01  1.8400659386536070e-01
3.0000000000000000e+00  1.3968082308703464e+00  8.4604444365799358e-01  1.6136031954087834e-01
3.2000000000000002e+00  1.5690949600067619e+00  8.7608145517247749e-01  1.3912805557242577e-01
3.4000000000000004e+00  1.7469500939488432e+00  9.0176122143676474e-01  1.1787624608752340e-01
3.6000000000000001e+00  1.9295251697605302e+00  9.2332966588164878e-01  9.8086278784280084e-02
3.8000000000000003e+00  2.1160298171720697e+00  9.4111799672931340e-01  8.0125918139197061e-02
4.0000000000000000e+00  2.3057464184620726e+00  9.5551822981069434e-01  6.4234121091690577e-02
4.2000000000000002e+00  2.4980396627069701e+00  9.6695707375360584e-01  5.0519747486645686e-02
4.4000000000000004e+00  2.6923609375430866e+00  9.7587083213647785e-01  3.8972610853629498e-02
4.6000000000000005e+00  2.8882479900178160e+00  9.8268350076070032e-01  2.9483772011648642e-02
4.8000000000000007e+00  3.0853206551774823e+00  9.8778952620078742e-01  2.1871186347443238e-02
5.0000000000000000e+00  3.2832736651563144e+00  9.9154190016439347e-01  1.5906798685318153e-02
5.2000000000000002e+00  3.4818676115086635e+00  9.9424553535992810e-01  1.1341788968929135e-02
5.4000000000000004e+00  3.6809190628975084e+00  9.9615530396275309e-01  7.9276598147061551e-03
5.6000000000000005e+00  3.8802906776333757e+00  9.9747776821291501e-01  5.4319576799275147e-03
5.8000000000000007e+00  4.0798819392396917e+00  9.9837549365019085e-01  3.6484136667473462e-03
6.0000000000000000e+00  4.2796209225138426e+00  9.9897287243586053e-01  2.4020398437572996e-03
6.2000000000000002e+00  4.4794572972826767e+00  9.9936254171909067e-01  1.5501706906563828e-03
6.4000000000000004e+00  4.6793566154314172e+00  9.9961170171767588e-01  9.8061511700977232e-04
6.6000000000000005e+00  4.8792958110602429e+00  9.9976787021003244e-01  6.0804426478449074e-04
6.8000000000000007e+00  5.0792597724490030e+00  9.9986381903703803e-01  3.6956257014031428e-04
7.0000000000000000e+00  5.2792388110290958e+00  9.9992160414794995e-01  2.2016895527113460e-04
7.2000000000000002e+00  5.4792268473431696e+00  9.9995571727792665e-01  1.2856980723515674e-04
7.4000000000000004e+00  5.6792201473338295e+00  9.9997545768489782e-01  7.3592983389223794e-05
7.6000000000000005e+00  5.8792164658040509e+00  9.9998665513912843e-01  4.1290311113698544e-05
7.8000000000000007e+00  6.0792144810719346e+00  9.9999288116610119e-01  2.2707751402803766e-05
8.0000000000000000e+00  6.2792134313460615e+00  9.9999627453530060e-01  1.2240926243249920e-05
8.2000000000000011e+00  6.4792128866784848e+00  9.9999808745923757e-01  6.4679786108453658e-06
8.4000000000000004e+00  6.6792126094414330e+00  9.9999903687484326e-01  3.3499397531974408e-06
8.5999999999999996e+00  6.8792124710150580e+00  9.9999952424780447e-01  1.7006679885717248e-06
8.8000000000000007e+00  7.0792124032167214e+00  9.9999976948975089e-01  8.4628412130585306e-07
9.0000000000000000e+00  7.2792123706452294e+00  9.9999989045379900e-01  4.1278790156475232e-07
9.2000000000000011e+00  7.4792123552968919e+00  9.9999994893901312e-01  1.9735668380412326e-07
9.4000000000000004e+00  7.6792123482031105e+00  9.9999997665715090e-01  9.2489158677487331e-08
9.6000000000000014e+00  7.8792123449874136e+00  9.9999998953401714e-01  4.2485812620375726e-08
9.8000000000000007e+00  8.0792123435577192e+00  9.9999999539788897e-01  1.9129831304430804e-08
1.0000000000000000e+01  8.2792123429343132e+00  9.9999999801539019e-01  8.4429158651022967e-09
1.0200000000000001e+01  8.4792123426677222e+00  9.9999999916068538e-01  3.6524803913092198e-09
1.0400000000000000e+01  8.6792123425559158e+00  9.9999999965190545e-01  1.5488074704694815e-09
1.0600000000000001e+01  8.8792123425099305e+00  9.9999999985842569e-01  6.4375569794427733e-10
1.0800000000000001e+01  9.0792123424913829e+00  9.9999999994353506e-01  2.6227617757128940e-10
1.1000000000000000e+01  9.2792123424840440e+00  9.9999999997791622e-01  1.0473955128094568e-10
1.1200000000000001e+01  9.4792123424811994e+00  9.9999999999153044e-01  4.0999322474265197e-11
1.1400000000000000e+01  9.6792123424801169e+00  9.9999999999681477e-01  1.5731015212516414e-11
1.1600000000000001e+01  9.8792123424797147e+00  9.9999999999882527e-01  5.9163099772618016e-12
1.1800000000000001e+01  1.0079212342479567e+01  9.9999999999957512e-01  2.1810177063481187e-12
1.2000000000000000e+01  1.0279212342479513e+01  9.9999999999984923e-01  7.8810042675859475e-13
1.2200000000000001e+01  1.0479212342479494e+01  9.9999999999994749e-01  2.7913740190928823e-13
1.2400000000000000e+01  1.0679212342479486e+01  9.9999999999998201e-01  9.6910000670083129e-14
1.2600000000000001e+01  1.0879212342479486e+01  9.9999999999999389e-01  3.2978678842745192e-14
1.2800000000000001e+01  1.1079212342479485e+01  9.9999999999999789e-01  1.1000489193601443e-14
1.3000000000000000e+01  1.1279212342479484e+01  9.9999999999999922e-01  3.5967050786090774e-15
1.3200000000000001e+01  1.1479212342479485e+01  9.9999999999999967e-01  1.1526879056075572e-15
1.3400000000000000e+01  1.1679212342479484e+01  9.9999999999999978e-01  3.6210349563133137e-16
1.3600000000000000e+01  1.1879212342479484e+01  9.9999999999999978e-01  1.1149817624085130e-16

If you look at the last two rows of fp you'll see that double precision has most likely run out of steam by eta = 13.4. Using gsl_odeiv2_step_rk4imp or some other choice instead of gsl_odeiv2_step_rk8pd will tend to cause the solver to fail somewhere before here rather than to simply stop producing updated values in fp. Interestingly, gsl_odeiv2_step_rk8pd will run out to about eta = 38.6 producing fpp = 1.35e-155 before fpp becomes non-decreasing (indicating the physics is glaringly b0rked as opposed to just wholly b0rked).

For want of a 2-line patch...

There's been much brouhaha about the Dread Pirate Roberts getting partially nabbed for using Stack Overflow. What struck me most about his question is he asked it only because PHP was missing a numeric constant, CURLPROXY_SOCKS5_HOSTNAME, which it should have had. Sure enough, there's an unaddressed bug (as of PHP 5.5.4) about the missing constant with an attached patch:

Developer: dr.scre@yandex.com

--- interface.c 2013-08-16 00:42:04.000000000 +0400
+++ interface_my.c 2013-08-18 03:12:28.948548811 +0400
@@ -827,6 +827,8 @@
  REGISTER_CURL_CONSTANT(CURLPROXY_HTTP);
  REGISTER_CURL_CONSTANT(CURLPROXY_SOCKS4);
  REGISTER_CURL_CONSTANT(CURLPROXY_SOCKS5);
+ REGISTER_CURL_CONSTANT(CURLPROXY_SOCKS4A);
+ REGISTER_CURL_CONSTANT(CURLPROXY_SOCKS5_HOSTNAME);
 
  /* Curl Share constants */
  REGISTER_CURL_CONSTANT(CURLSHOPT_NONE);

For want of a 2-line patch, a $1.2B empire was lost.

On Working in the Dark.

At 7 PM tonight, despite my father's long-standing advice about working on things when pressed for time, I start futzing with the right turn signal switch on a recently-new-to-me '99 BMW R1100S. The switch has been skittish the last two days. It's no fun to find out that traffic doesn't think you're switching lanes when you think you've been clearly signalling it for an eighth of a mile. Michelle has to leave to teach yoga at 7:45 PM. This leaves me forty-five minutes to work before I must watch the kiddo. Night is about to fall. Sans garage, I'm working in the driveway.

I have the handlebar control assembly down to the switch innards apart in five minutes, see some grime, and decide to make a run for some electrical cleaner rather than to try my luck applying to the contact surfaces the fine-grit sandpaper that I brought outside. To give Michelle a chance to get ready to teach, I snag the two-and-a-half-year-old in his PJs and we go by car to an auto store on East 7th. After a long, slow line, we're back at 7:30. I get the switch cleaned out (after accidentally shooting electrical cleaner into a hangnail. Don't. Ever. Do. This.) and start putting stuff back together.

Now it's getting dark. Michelle brings out a flashlight. I get the electrical bits back together. Michelle has to leave and brings out our shoeless son Ozark, his Magna Doodle, and his wooden toy pliers. She coaxes Ozark into staying put next to me on the driveway and takes off. He miraculously stays out of the small collection of bike parts and is content to play around on the driveway while I chew on a flashlight trying to point it in the right direction as I re-insert screws and curse the finicky-to-button-up assembly. Ozark disappears over to the left side of the bike and is farting around with his wooden pliers near the left fork. I ask him to come around to my side and he ignores me. No matter as he's being uncommonly good. I periodically inquire what he's up to and he responds promptly.

At 7:55 I get everything back together well-enough that the repair will sit for the night and I can double check it in the morning. I'm elated that nothing bad happened on a hurry-up repair job in the dusk. This has never happened to me in my life.

Suddenly, I hear a weird scritching sound. I ask Ozark what he's doing. Out of the dark on the far side of the bike a little voice proudly replies "I'm sanding the motorcycle". He's got the fine-grit sandpaper and is rubbing it against the fairing near the front headlight. Fortunately the cosmetic damage is light. I crack up. He cracks up. I call my dad to convey the story. Ozark re-cracks up when he hears my dad crack up on the phone. Then the kiddo and I go inside and drink far too much chocolate milk.

An iostream insertion wrapper to provide full-precision floating point output for C++

In C++ one often wants output full precision floating point data. Iostreams make this a pain. Especially when you'd like the results to line up nicely in your terminal without requiring everything to be in scientific notation. Because reading "12" as 12 is easier than parsing 1.2e+01.

From my research code Suzerain, here's a header-only stream insertion wrapper up to the task. First, a simple use case with output:

using suzerain::fullprec;
cout << fullprec<>(1.23456789012345678e-5) << endl; // "   1.23456789012e-05"
cout << fullprec<>(1.23456789012345678e-4) << endl; // "   1.23456789012e-04"
cout << fullprec<>(1.23456789012345678e-3) << endl; // "   1.23456789012e-03"
cout << fullprec<>(1.23456789012345678e-2) << endl; // "   0.012345678901235"
cout << fullprec<>(1.23456789012345678e-1) << endl; // "   0.123456789012346"
cout << fullprec<>(1.23456789012345678e+0) << endl; // "   1.234567890123457"
cout << fullprec<>(1.23456789012345678e+1) << endl; // "  12.345678901234567"
cout << fullprec<>(1.23456789012345678e+2) << endl; // " 123.456789012345681"
cout << fullprec<>(1.23456789012345678e+3) << endl; // "   1.23456789012e+03"
cout << fullprec<>(1.23456789012345678e+4) << endl; // "   1.23456789012e+04"
cout << fullprec<>(1.23456789012345678e+5) << endl; // "   1.23456789012e+05"

Now the code. You'll need to swap real_t for double in your context as suzerain::real_t is defined elsewhere.

Worst open source release ever: Suzerain

Though it's far from in a public-ready state, I'm tired of my research code, Suzerain, having no public presence whatsoever:

• Wiki
• Roadmap
• Issue Tracking
• Bitbucket Repository (though a word of caution-- the main repository is still private and this is a mirror)
• GitHub Repository (ditto, because who can pick between 'em, right?)
• Google Group

Though primarily my baby up until a year ago, the cast of characters now includes all of whom are affiliated with The Center for Predictive Engineering and Computational Science (PECOS) within The Institute for Computational Engineering and Sciences (ICES) at The University of Texas at Austin.

This work is supported by the Department of Energy [National Nuclear Security Administration] under Award Number [DE-FC52-08NA28615].

Somewhat taming the bump function's derivatives

Often I need to initialize some smooth pulse in a simulation. Today, for example, I'm working with acoustic and entropy pulses to test the effectiveness of nonreflecting boundary conditions for the Euler equations (see, e.g., Baum et al). Bump functions are nice for this purpose as they're infinitely smooth and have compact support.

Awhile back I realized that using just a classical bump function produced large gradients which had to be resolved. Years ago I used Mathematica to see if taking some power of the classical bump function, which remains infinitely smooth, could produce smaller gradients better suited for numerical use.

The punchline is that using the fourth power of the bump function, BumpPower[4] in the following parlance, is Goodness (TM). The following Mathematica notebook demonstrates the problem and discusses how one arrives at this empirical result:

The obligatory XMonad configuration post

Since I've been loving XMonad for the past several months on a 2-up desktop and a 1-up laptop, here's the obligatory personal XMonad configuration post:

General comments:

1. This started from the vanilla, session-manager-aware Fedora XMonad configuration.
2. The Windows key is the modifier.
3. The default dmenu_run key, mod-p, interferes with some monitor-swapping hotkey on Ubuntu. Hence the second keymapping.
4. Toggling maximization in one keypress has been handy.
5. I like the mod+arrow combinations better than the suggested CycleWS ones. Rather than moving through a bunch of empty workspaces, these permit you to just move through workspaces with windows and to get a clean workspace with a single keypress. Much thanks to Marshall Lochbaum and Brandon Allbery for chiming in on the XMonad mailing list with the followTo implementation.
6. The prev/next screen items permit throwing windows back and forth when 2-up.
7. I find the REFLECTX binding handy when 2-up to push the non-main chat window off to the leftmost monitor position.
8. Layout.Tabbed's simpleTabbed is additionally available as a layout option.
9. New windows spawn outside the master pane with the master window keeping focus.

Please drop a comment if you happen to find anything useful or want to contribute a tweak.

Python should have Haskell's $ operator

In Python, I often want to repeatedly say "apply some function to the remainder of the current statement". Think foo(bar(baz(qux, quux))). In my text editor, which mostly automatically handles the parenthesis correctly, this is a quick thing to type.

Within IPython, statements like that are tedious. IPython does have a limited version of this in its automatic parenthesis feature. Notice how many syntax error gotchas lurk in those examples.

Haskell's $ operator is pretty sweet. All of the following are equivalent in Haskell:

putStrLn (show (1 + 1))
putStrLn (show $ 1 + 1)
putStrLn $ show (1 + 1)
putStrLn $ show $ 1 + 1

Python should adopt Haskell's $ operator. That is, the following should be equivalent:

foo(bar(baz(qux, quux)))
foo $ bar $ baz(qux, quux)
foo $ bar $ baz $ qux, quux

Though odd looking, there's easy consistency for unary functions:

foo()
foo $

Haters of the new-ish print might find respite in print $ foo.

Plus ça change.org

This past fall I spun up a change.org petition entitled 83rd Texas State Legislature: Reduce traffic congestion by legalizing lane-splitting for motorcycles. As of right now, 1,664 people signed the petition. Sadly, despite gathering this tome of names and much late-night letter writing to both newspapers and elected officials, the Texas legislature didn't introduce the desired bill.

But there may be hope for the 84th legislature in two years... Today I got an unexpected and very cool email from the office of Senator Kirk Watson. Watson said he'll try to get lane splitting formally studied between now and the next legislature. I'll redefine success to match this outcome and claim it. The full message follows.

How to not land a job.

Quick benchmark on the bounded_buffer<T> example from Boost Templated Circular Buffer

Boost's Templated Circular Buffer documentation claims that some modest Boost.Thread wrapping of a boost::circular_buffer<T> makes a bounded producer-consumer queue. Happily, they provide a benchmark. I believe 'em, of course, but I'd like to get an idea of the variability behind the claims.

Rusty Russell, of CCAN fame, put out a nice pipe-oriented processor for gathering statistics. This can directly be fed the output of the benchmark as a pretty cute one-liner:

gcc -Wall -O3 -pthread bounded_buffer_comparison.cpp -lboost_thread-mt
for i in `seq 1 100`; do ./a.out; done | ./stats

With a touch of reformatting and sorting the columns by increasing mean, here are the results on a circa-2007 T61 stinkpad of mine:

bounded_buffer<int>                         0.45 -- 0.98  (0.7588 +/- 0.17 ) s
bounded_buffer<std::string>                 0.58 -- 0.91  (0.7678 +/- 0.095) s
bounded_buffer_deque_based<std::string>     0.43 -- 0.94  (0.8291 +/- 0.095) s
bounded_buffer_deque_based<int>             0.48 -- 1.03  (0.8374 +/- 0.12 ) s
bounded_buffer_space_optimized<int>         0.55 -- 1.10  (0.8459 +/- 0.1  ) s
bounded_buffer_space_optimized<std::string> 0.46 -- 1.21  (0.8852 +/- 0.24 ) s
bounded_buffer_list_based<int>              2.84 -- 6.65  (4.2641 +/- 0.74 ) s
bounded_buffer_list_based<std::string>      2.22 -- 5.77  (4.2663 +/- 0.66 ) s

Building -O3 makes a big difference as you might expected and I've not bothered showing -O2. On average, the sample bounded_buffer<T> does slightly outperform the deque-based logic as the Boost Templated Circular Buffer documentation suggests. The variability is, uhhh, more variable than I'd expected.

Mergeable accumulation of the running min, mean, max, and variance

Boost Accumulators are pretty slick, but they do have the occasional shortcoming. To my knowledge, one cannot merge data from independent accumulator sets. Merging is a nice operation to have when you want to compute statistics in parallel and then roll up the results into a single summary.

In his article Accurately computing running variance, John D. Cook presents a small class for computing a running mean and variance using an algorithm with favorable numerical properties reported by Knuth in TAOCP. Departing from Cook's sample, my rewrite below includes

• templating on the floating point type,
• permitting accumulating multiple statistics simultaneously without requiring multiple counters,
• additionally tracking the minimum and maximum while holding storage overhead constant relative to Cook's version,
• providing sane (i.e. NaN) behavior when no data has been processed,
• permitting merging information from multiple instances,
• permitting clearing an instance for re-use,
• assertions to catch when one screws up,
• and adding Doxygen-based documentation.

Please let me know if you find this useful. If anyone is interested, I will also post the associated regression tests.

bash parameter substitution quick reference

I can never remember the syntax for all the various substitution possibilities in bash. Though many places detail them all (e.g. man bash, Advanced Bash Scripting Guide), the syntax is threaded through paragraphs of text. Here is my distilled quick reference.

Many thanks to snipt for housing this and the many other miscellaneous code snippets I can never keep in my head but constantly need.

Candidacy!

This past Monday I finally made candidacy in my graduate program.