Avoiding tedious numerics using Boost Accumulators

For a turbulence warm up assignment, we needed to compute some correlation coefficients for trajectory behavior on the Lorenz attractor. Specifically, compute the mean, variance, and cross-correlations of the three solution components. I'd skimmed through the Boost Accumulator Framework a couple of months back and thought it looked interesting:

#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/covariance.hpp>
#include <boost/accumulators/statistics/mean.hpp>
#include <boost/accumulators/statistics/stats.hpp>
#include <boost/accumulators/statistics/variance.hpp>
#include <boost/accumulators/statistics/variates/covariate.hpp>
#include <boost/format.hpp>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>

// Read an input streams and compute the relevant statistics
template<typename charT, typename traits>
int process(std::basic_ostream<charT, traits>& out,
            std::basic_istream<charT, traits>& in) {

    using namespace boost::accumulators;
    accumulator_set<double, stats<tag::mean,
                                  tag::variance,
                                  tag::covariance<double, tag::covariate1> > > x_acc;
    accumulator_set<double, stats<tag::mean,
                                  tag::variance,
                                  tag::covariance<double, tag::covariate1> > > y_acc;
    accumulator_set<double, stats<tag::mean,
                                  tag::variance,
                                  tag::covariance<double, tag::covariate1> > > z_acc;

    double t, x, y, z;
    std::basic_string<charT, traits> line;
    boost::format fmt(" %016e   %016e   %016e   %016e   %016e   %016e   %016e   %016e   %016e   %016e");
    while (in) {
        if (!std::getline(in, line)) break; // Read columns per line to
        std::istringstream in(line);       //    be a little defensive
        in >> t >> x >> y >> z;

        // Compute running statistics
        x_acc(x, covariate1 = y);
        y_acc(y, covariate1 = z);
        z_acc(z, covariate1 = x);

        // Output running statistics
        out << fmt % t
                   % mean(x_acc)
                   % mean(y_acc)
                   % mean(z_acc)
                   % variance(x_acc)
                   % variance(y_acc)
                   % variance(z_acc)
                   % covariance(x_acc)
                   % covariance(z_acc)
                   % covariance(y_acc)
            << std::endl;
    }

    return 0;
}

Compared with the alternative loops/bookkeeping/etc., this is a pretty slick alternative. The only downside was that the accumulator declarations took me awhile to decipher from the documentation.

As you probably noticed, I am reading/writing plain text data files and computing the statistics at every time step. This is overkill and slow. I'd probably have fixed that by now if it weren't for the fact that I ran across Pipe Viewer in the last two days. Quite a neat shell tool that.