ar.hpp

Go to the documentation of this file.

// Copyright (C) 2012, 2013, 2026 Rhys Ulerich
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
#ifndef AR_HPP
#define AR_HPP
 
#ifndef AR_SUPPRESS_DOXYGEN_MAINPAGE
 
#endif /* AR_SUPPRESS_DOXYGEN_MAINPAGE */
 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <functional>
#include <iterator>
#include <limits>
#include <numeric>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace ar
{
 
 
#define AR_GCC_VERSION \
    (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
 
#define AR_STRINGIFY_HELPER(x) #x
 
#define AR_STRINGIFY(x) AR_STRINGIFY_HELPER(x)
 
#define AR_ENSURE_MSGEXCEPT(expr, msg, except) \
    if (!(expr)) throw except(msg)
 
#define AR_ENSURE_MSG(expr, msg) \
    AR_ENSURE_MSGEXCEPT(expr, msg, std::logic_error)
 
#define AR_ENSURE(expr) \
    AR_ENSURE_MSG(expr, AR_STRINGIFY(expr)" false")
 
#define AR_ENSURE_ARG(expr) \
    AR_ENSURE_MSGEXCEPT(expr, AR_STRINGIFY(expr)" false", std::invalid_argument)
 
#define AR_ENSURE_EXCEPT(expr, except) \
    AR_ENSURE_MSGEXCEPT(expr, AR_STRINGIFY(expr)" false", except)
 
 
template <typename InputIterator,
          typename OutputType1,
          typename OutputType2>
std::size_t welford_nvariance(InputIterator first,
                              InputIterator last,
                              OutputType1&   mean,
                              OutputType2&   nvar)
{
    using std::iterator_traits;
    using std::size_t;
    typedef typename iterator_traits<InputIterator>::value_type value;
 
    size_t N  = 1;  // Running next sample number
    value  m  = 0;  // Running mean of data thus far
    value  nv = 0;  // Running variance times the number of samples
 
    while (first != last)
    {
        value x  = *first++;
        value d  = x - m;
        m       += d / N++;
        nv      += d*(x - m);
    }
 
    mean = m;
    nvar = nv;
    return N-1;
}
 
template <typename InputIterator,
          typename OutputType1,
          typename OutputType2>
std::size_t welford_variance_population(InputIterator first,
                                        InputIterator last,
                                        OutputType1&  mean,
                                        OutputType2&  var)
{
    using std::size_t;
    size_t N = welford_nvariance(first, last, mean, var);
    var /= N;
    return N;
}
 
template <typename InputIterator,
          typename OutputType1,
          typename OutputType2>
std::size_t welford_variance_sample(InputIterator first,
                                    InputIterator last,
                                    OutputType1&  mean,
                                    OutputType2&  var)
{
    using std::size_t;
    size_t N = welford_nvariance(first, last, mean, var);
    var /= (N - 1);
    return N;
}
 
template <typename InputIterator1,
          typename InputIterator2,
          typename OutputType1,
          typename OutputType2,
          typename OutputType3>
std::size_t welford_ncovariance(InputIterator1 first1,
                                InputIterator1 last1,
                                InputIterator2 first2,
                                OutputType1&   mean1,
                                OutputType2&   mean2,
                                OutputType3&   ncovar)
{
    using std::iterator_traits;
    using std::size_t;
    typedef typename iterator_traits<InputIterator1>::value_type value1;
    typedef typename iterator_traits<InputIterator2>::value_type value2;
 
    size_t      N  = 1;  // Running next sample number
    value1      m1 = 0;  // Running mean of first data set thus far
    value2      m2 = 0;  // Running mean of second data set thus far
    OutputType3 nc = 0;  // Running covariance times the number of samples
 
    while (first1 != last1)
    {
        value1 x1  = *first1++;
        value1 d1  = x1 - m1;
        m1        += d1 / N;
 
        value2 x2  = *first2++;
        value2 d2  = x2 - m2;
        m2        += d2 / N;
 
        nc += d1*(x2 - m2);
 
        ++N;
    }
 
    mean1  = m1;
    mean2  = m2;
    ncovar = nc;
    return N-1;
}
 
template <typename InputIterator1,
          typename InputIterator2,
          typename OutputType1,
          typename OutputType2,
          typename OutputType3>
std::size_t welford_covariance_population(InputIterator1 first1,
                                          InputIterator1 last1,
                                          InputIterator2 first2,
                                          OutputType1&   mean1,
                                          OutputType2&   mean2,
                                          OutputType3&   covar)
{
    using std::size_t;
    size_t N = welford_ncovariance(first1, last1, first2, mean1, mean2, covar);
    covar /= N;
    return N;
}
 
template <typename InputIterator1,
          typename InputIterator2,
          typename OutputType1,
          typename OutputType2,
          typename OutputType3>
std::size_t welford_covariance_sample(InputIterator1 first1,
                                      InputIterator1 last1,
                                      InputIterator2 first2,
                                      OutputType1&   mean1,
                                      OutputType2&   mean2,
                                      OutputType3&   covar)
{
    using std::size_t;
    size_t N = welford_ncovariance(first1, last1, first2, mean1, mean2, covar);
    covar /= (N - 1);
    return N;
}
 
template <typename InputIterator,
          typename ValueType>
ValueType welford_inner_product(InputIterator first,
                                InputIterator last,
                                ValueType     init)
{
    typename std::iterator_traits<InputIterator>::value_type mean;
    ValueType nvar;
    const std::size_t N = welford_nvariance(first, last, mean, nvar);
    return init + (nvar + N*(mean*mean));
}
 
template <typename InputIterator1,
          typename InputIterator2,
          typename ValueType>
ValueType welford_inner_product(InputIterator1 first1,
                                InputIterator1 last1,
                                InputIterator2 first2,
                                ValueType      init)
{
    typename std::iterator_traits<InputIterator1>::value_type mean1;
    typename std::iterator_traits<InputIterator2>::value_type mean2;
    ValueType ncovar;
    const std::size_t N = welford_ncovariance(
            first1, last1, first2, mean1, mean2, ncovar);
    return init + (ncovar + N*(mean1*mean2));
}
 
 
#if _MSC_VER > 1400
# pragma float_control(push)
# pragma float_control(precise, on)
#endif
 
template <typename ValueType,
          typename InputIterator1,
          typename InputIterator2>
ValueType
#if (AR_GCC_VERSION > 40305)
    __attribute__((__optimize__("no-associative-math")))
#endif
negative_half_reflection_coefficient(InputIterator1 a_first,
                                     InputIterator1 a_last,
                                     InputIterator2 b_first)
#if (AR_GCC_VERSION > 40305) ||  (_MSC_VER > 1400)
{
    ValueType ns = 0, nt, nc = 0, ny;  // Kahan numerator accumulation
    ValueType ds = 0, dt, dc = 0, dy;  // Kahan denominator accumulation
 
    while (a_first != a_last)
    {
        ValueType xa = *a_first++;     // Denominator: \vec{a}\cdot\vec{a}
        dy = (xa * xa) - dc;
        dt = ds + dy;
        dc = (dt - ds) - dy;
        ds = dt;
 
        ValueType xb = *b_first++;     // Denominator: \vec{b}\cdot\vec{b}
        dy = (xb * xb) - dc;
        dt = ds + dy;
        dc = (dt - ds) - dy;
        ds = dt;
 
        ny = (xa * xb) - nc;           // Numerator:   \vec{a}\cdot\vec{b}
        nt = ns + ny;
        nc = (nt - ns) - ny;
        ns = nt;
    }
 
    return ns + nc == 0                // Does special zero case apply?
        ? 0                            // Yes, to avoid NaN from 0 / 0
        : (ns + nc) / (ds + dc);       // No, correct final sums and form ratio
}
#else
#warning Using Non-Kahan version of ar::negative_half_reflection_coefficient.
{
    ValueType ns = 0;
    ValueType ds = 0;
 
    while (a_first != a_last)
    {
        ValueType xa  = *a_first++;
        ValueType xb  = *b_first++;
        ns           += xa * xb;            // Numerator
        ds           += xa * xa + xb * xb;  // Denominator
    }
 
    return ns == 0                          // Does special zero case apply?
        ? 0                                 // Yes, to avoid NaN from 0 / 0
        : ns / ds;                          // No, form ratio
}
#endif
 
#if _MSC_VER > 1400
# pragma float_control(pop)
#endif
 
template <class InputIterator,
          class Value,
          class OutputIterator1,
          class OutputIterator2,
          class OutputIterator3,
          class OutputIterator4,
          class Vector>
std::size_t burg_method(InputIterator   data_first,
                        InputIterator   data_last,
                        Value&          mean,
                        std::size_t&    maxorder,
                        OutputIterator1 params_first,
                        OutputIterator2 sigma2e_first,
                        OutputIterator3 gain_first,
                        OutputIterator4 autocor_first,
                        const bool      subtract_mean,
                        const bool      hierarchy,
                        Vector&         f,
                        Vector&         b,
                        Vector&         Ak,
                        Vector&         ac)
{
    using std::bind2nd;
    using std::copy;
    using std::distance;
    using std::fill;
    using std::inner_product;
    using std::min;
    using std::minus;
    using std::size_t;
 
    // Initialize f from [data_first, data_last) and fix number of samples
    f.assign(data_first, data_last);
    const size_t N = f.size();
 
    // Stably compute the incoming data's mean and population variance
    mean = 0;
    Value sigma2e = 0;
    welford_variance_population(f.begin(), f.end(), mean, sigma2e);
 
    // When requested, subtract the just-computed mean from the data.
    // Adjust, if necessary, to make sigma2e the second moment.
    if (subtract_mean)
    {
        for (typename Vector::iterator it = f.begin(); it != f.end(); ++it) {
            *it -= mean;
        }
    }
    else
    {
        sigma2e += mean*mean;
    }
 
    // At most maxorder N-1 can be fit from N samples.  Beware N is unsigned.
    maxorder = (N == 0) ? 0 : min(static_cast<size_t>(maxorder), N-1);
 
    // Output sigma2e and gain for a zeroth order model, if requested.
    Value gain = 1;
    if (hierarchy || maxorder == 0)
    {
        *sigma2e_first++ = sigma2e;
        *gain_first++    = gain;
    }
 
    // Initialize and perform Burg recursion
    if (maxorder) b = f;  // Copy iff non-trivial work required
    Ak.assign(maxorder + 1, Value(0));
    Ak[0] = 1;
    ac.clear();
    ac.reserve(maxorder);
    for (size_t kp1 = 1; kp1 <= maxorder; ++kp1)
    {
        // Compute mu from f, b, and Dk and then update sigma2e and Ak using mu
        // Afterwards, Ak[1:kp1] contains AR(k) coefficients by the recurrence
        // Must treat mu result of 0 / 0 as 0 to avoid NaNs on constant signals
        // By the recurrence, Ak[kp1] will also be the reflection coefficient
        Value mu = -2 * negative_half_reflection_coefficient<Value>(
                f.begin() + kp1, f.end(), b.begin());
 
        sigma2e *= (1 - mu*mu);
        for (size_t n = 0; n <= kp1/2; ++n)
        {
            Value t1 = Ak[n] + mu*Ak[kp1 - n];
            Value t2 = Ak[kp1 - n] + mu*Ak[n];
            Ak[n] = t1;
            Ak[kp1 - n] = t2;
        }
 
        // Update the gain per Broersen 2006 equation (5.25)
        gain *= 1 / (1 - Ak[kp1]*Ak[kp1]);
 
        // Compute and output the next autocorrelation coefficient
        // See Broersen 2006 equations (5.28) and (5.31) for details
        ac.push_back(-inner_product(ac.rbegin(), ac.rend(),
                                    Ak.begin() + 1, Ak[kp1]));
 
        // Output parameters and the input and output variances when requested
        if (hierarchy || kp1 == maxorder)
        {
            params_first = copy(Ak.begin() + 1, Ak.begin() + kp1 + 1,
                                params_first);
            *sigma2e_first++ = sigma2e;
            *gain_first++    = gain;
        }
 
        // Update f and b for the next iteration if another remains
        if (kp1 < maxorder)
        {
            for (size_t n = 0; n < N - kp1; ++n)
            {
                Value t1 = f[n + kp1] + mu*b[n];
                Value t2 = b[n] + mu*f[n + kp1];
                f[n + kp1] = t1;
                b[n] = t2;
            }
        }
    }
 
    // Output the lag [0,maxorder] autocorrelation coefficients in single pass
    *autocor_first++ = 1;
    copy(ac.begin(), ac.end(), autocor_first);
 
    // Return the number of values processed in [data_first, data_last)
    return N;
}
 
template <class InputIterator,
          class Value,
          class OutputIterator1,
          class OutputIterator2,
          class OutputIterator3,
          class OutputIterator4>
std::size_t burg_method(InputIterator     data_first,
                        InputIterator     data_last,
                        Value&            mean,
                        std::size_t&      maxorder,
                        OutputIterator1   params_first,
                        OutputIterator2   sigma2e_first,
                        OutputIterator3   gain_first,
                        OutputIterator4   autocor_first,
                        const bool        subtract_mean = false,
                        const bool        hierarchy     = false)
{
    using std::vector;
    vector<Value> f, b, Ak, ac; // Working storage
 
    return burg_method(data_first, data_last, mean, maxorder,
                       params_first, sigma2e_first, gain_first,
                       autocor_first, subtract_mean, hierarchy,
                       f, b, Ak, ac);
}
 
// Type erasure for NoiseGenerator parameters within predictor.
// Either std::tr1::function or boost::function would better provide the
// desired capability but both add additional, undesired dependencies.
namespace
{
 
template <typename Value>
struct nullary
{
    virtual ~nullary() {}
    virtual Value operator()() = 0;
    virtual nullary* clone()   = 0;
};
 
template<typename Value>
struct nullary_impl0 : public nullary<Value>
{
    Value operator()()
    {
        return 0;
    }
    nullary_impl0* clone()
    {
        return new nullary_impl0();
    }
};
 
template<typename Value, class T>
struct nullary_impl1 : public nullary<Value>
{
    nullary_impl1(T t) : t(t) {}
    Value operator()()
    {
        return t();
    }
    nullary_impl1* clone()
    {
        return new nullary_impl1(t);
    }
    T t;
};
 
}
 
template <typename Value, typename Index = std::size_t>
class predictor
{
public:
    typedef std::input_iterator_tag iterator_category;
    typedef Value                   value_type;
    typedef std::ptrdiff_t          difference_type;
    typedef const Value*            pointer;
    typedef const Value&            reference;
 
    explicit predictor(Index n = 0) : n(n), d(), g(0), xn()
    {
#ifndef NDEBUG
        using std::numeric_limits;
        if (numeric_limits<Value>::has_quiet_NaN)
            xn = numeric_limits<Value>::quiet_NaN();
#endif
    }
 
    template <class RandomAccessIterator>
    predictor(RandomAccessIterator params_first,
              RandomAccessIterator params_last)
        : n(0),
          d(2*std::distance(params_first, params_last), 0),
          g(new nullary_impl0<Value>()),
          xn((*g)())
    {
        // Finish preparing d = [ a_p, ..., a_1, 0, ..., 0 ]
        using std::vector;
        typename vector<Value>::size_type i = d.size() / 2;
        while (i --> 0) d[i] = *params_first++;
 
        // Now x_n = 0 because x_{n-p} = ... = x_{n-1} = 0 by construction.
    }
 
    template <class RandomAccessIterator,
              class NoiseGenerator>
    predictor(RandomAccessIterator params_first,
              RandomAccessIterator params_last,
              NoiseGenerator generator)
        : n(0),
          d(2*std::distance(params_first, params_last), 0),
          g(new nullary_impl1<Value,NoiseGenerator>(generator)),
          xn((*g)())
    {
        // Finish preparing d = [ a_p, ..., a_1, 0, ..., 0 ]
        using std::vector;
        typename vector<Value>::size_type i = d.size() / 2;
        while (i --> 0) d[i] = *params_first++;
 
        // Here x_0 = \epsilon_0 because x_{0-p} = ... = x_{0-1} = 0.
    }
 
    predictor(const predictor& other)
        : n(other.n),
          d(other.d),
          g(other.g ? other.g->clone() : 0),
          xn(other.xn)
    {}
 
    predictor& operator=(const predictor& other)
    {
        if (this != &other)
        {
            nullary<Value> *tmp = 0;
            try
            {
                tmp = other.g ? other.g->clone() : 0;
            }
            catch (...)
            {
                delete tmp;
                throw;
            }
            n = other.n;
            d = other.d;
            delete g;
            g = tmp;
            xn = other.xn;
        }
        return *this;
    }
 
    ~predictor()
    {
        delete g;
    }
 
    template <class InputIterator>
    predictor& initial_conditions(InputIterator initial_first,
                                  const Value x0adjust = 0)
    {
        // Zero the simulation time.
        n = 0;
 
        // Set d = [ a_p, ..., a_1, x_{n-p}, ..., x_{n-1} ]
        using std::vector;
        typename vector<Value>::size_type i = d.size();
        typename vector<Value>::size_type p = i / 2;
        while (i --> p) d[i] = *initial_first++;
 
        // Make x_n := - a_p*x_{n-p} - ... - a_1*x_{n-1} + x_n + x0adjust.
        // By design, x_n was whatever it happened to be.
        using std::inner_product;
        xn += x0adjust;
        xn  = -inner_product(d.begin(), d.begin() + p, d.begin() + p, -xn);
 
        return *this;
    }
 
    // Concept: InputIterator
 
    predictor& operator++()
    {
        using std::distance;
        using std::inner_product;
        using std::vector;
 
        if (g)
        {
            typename vector<Value>::size_type p = d.size() / 2;
            if (p)
            {
                // Make x_n = - a_p*x_{n-p} - ... - a_1*x_{n-1} + \epsilon_n
                // by (conceptually) storing previously computed x_n into
                // circular buffer, updating ++n, and computing x_{n+1}.
                typename vector<Value>::iterator ab = d.begin();
                typename vector<Value>::iterator xb = ab + p;
                typename vector<Value>::iterator c  = xb + n % p;
                typename vector<Value>::iterator xe = d.end();
                *c++ =  xn;
                xn   =  inner_product(c,  xe, ab,                   -(*g)());
                xn   = -inner_product(xb,  c, ab + distance(c, xe),   xn   );
            }
            else
            {
                xn = (*g)();
            }
        }
        else
        {
#ifndef NDEBUG
            using std::numeric_limits;
            if (numeric_limits<Value>::has_quiet_NaN)
                xn = numeric_limits<Value>::quiet_NaN();
#endif
        }
 
        ++n;
 
        return *this;
    }
 
    predictor operator++(int) // Postfix increment
    {
        predictor t(*this);
        ++*this;
        return t;
    }
 
    reference operator* () const
    {
        return xn;
    }
 
    // Concept: EqualityComparable
 
    bool operator== (const predictor& other) const
    {
        return n == other.n;
    }
 
    bool operator!= (const predictor& other) const
    {
        return !(*this == other);
    }
 
private:
 
    Index n;
 
    typename std::vector<Value> d;
 
    nullary<Value> *g;
 
    Value xn;
};
 
template <class RandomAccessIterator,
          class InputIterator,
          class Value>
predictor<typename std::iterator_traits<RandomAccessIterator>::value_type>
autocorrelation(RandomAccessIterator params_first,
                RandomAccessIterator params_last,
                Value                gain,
                InputIterator        autocor_first)
{
    predictor<typename std::iterator_traits<
            RandomAccessIterator
        >::value_type> p(params_first, params_last);
    p.initial_conditions(++autocor_first, 1 / gain);
    return p;
}
 
template <class Value>
Value decorrelation_time(const std::size_t N,
                         predictor<Value> rho,
                         const bool abs_rho = false)
{
    using std::abs;
    using std::size_t;
 
    Value T0 = *rho++;
 
    const Value twoinvN = Value(2) / N;
    if (abs_rho) {
        for (size_t i = 1; i <= N; ++i, ++rho)
            T0 += (2 - i*twoinvN) * abs(*rho);
    } else {
        for (size_t i = 1; i <= N; ++i, ++rho)
            T0 += (2 - i*twoinvN) * (*rho);
    }
 
    return T0;
}
 
template <class Value>
Value decorrelation_time(const std::size_t N,
                         predictor<Value> rho1,
                         predictor<Value> rho2,
                         const bool abs_rho = false)
{
    using std::abs;
    using std::size_t;
 
    Value T0 = 1;
    ++rho1;
    ++rho2;
 
    const Value twoinvN = Value(2) / N;
    if (abs_rho) {
        for (size_t i = 1; i <= N; ++i, ++rho1, ++rho2)
            T0 += (2 - i*twoinvN) * abs((*rho1) * (*rho2));
    } else {
        for (size_t i = 1; i <= N; ++i, ++rho1, ++rho2)
            T0 += (2 - i*twoinvN) * ((*rho1) * (*rho2));
    }
 
    return T0;
}
 
template<class RandomAccessIterator,
         class InputIterator,
         class OutputIterator>
void zohar_linear_solve(RandomAccessIterator a_first,
                        RandomAccessIterator a_last,
                        RandomAccessIterator r_first,
                        InputIterator        d_first,
                        OutputIterator       s_first)
{
    using std::copy;
    using std::distance;
    using std::inner_product;
    using std::invalid_argument;
    using std::iterator_traits;
    using std::reverse_iterator;
    using std::vector;
 
    // Tildes indicate transposes while hats indicate reversed vectors.
 
    // InputIterator::value_type determines the working precision
    typedef typename iterator_traits<InputIterator>::value_type value_type;
    typedef vector<value_type> vector_type;
    typedef typename vector_type::size_type size_type;
 
    // Determine problem size using [a_first,a_last) and ensure nontrivial
    typename iterator_traits<RandomAccessIterator>::difference_type dist
            = distance(a_first, a_last);
    if (dist < 1) throw invalid_argument("distance(a_first, a_last) < 1");
    const size_type n = static_cast<size_type>(dist);
 
    // Allocate working storage and set initial values for recursion:
    vector_type s;    s   .reserve(n+1); s   .push_back( *d_first);
    vector_type ehat; ehat.reserve(n  ); ehat.push_back(-a_first[0]);
    vector_type g;    g   .reserve(n  ); g   .push_back(-r_first[0]);
    value_type lambda  = 1 - a_first[0]*r_first[0];
 
    // Though recursive updates to s and g can be done in-place, updates to
    // ehat seemingly require one additional vector for storage:
    //
    // "This sequence of computations is economical of storage.  It is only
    // necessary to retain quantities computed at level m - 1 until the
    // computations at level m are complete." [Trench1967, page 1504]
    vector_type next_ehat; next_ehat.reserve(n);
 
    // Recursion for i = {1, 2, ..., n - 1}:
    for (size_type i = 1; i < n; ++i)
    {
 
        reverse_iterator<RandomAccessIterator> rhat_first(r_first + i);
 
        // \theta_i =  \delta_{i+1}  - \tilde{s}_i \hat{r}_i
        const value_type neg_theta = inner_product(
                s.begin(), s.end(), rhat_first, value_type(-(*++d_first)));
 
        // \eta_i   = -\rho_{-(i+1)} - \tilde{a}_i \hat{e}_i
        const value_type neg_eta   = inner_product(
                ehat.begin(), ehat.end(), a_first, value_type(a_first[i]));
 
        // \gamma_i = -\rho_{i+1}    - \tilde{g}_i \hat{r}_i
        const value_type neg_gamma = inner_product(
                g.begin(), g.end(), rhat_first, value_type(r_first[i]));
 
        /*
         * s_{i+1} = \bigl(\begin{smallmatrix}
         *              s_i + (\theta_i/\lambda_i) \hat{e}_i \\
         *              \theta_i/\lambda_i
         *          \end{smallmatrix}\bigr)
         *
         * \hat{e}_{i+1} = \bigl(\begin{smallmatrix}
         *                     \eta_i/\lambda_i \\
         *                     \hat{e}_i + (\eta_i/\lambda_i) g_i
         *                 \end{smallmatrix}\bigr)
         *
         * g_{i+1} = \bigl(\begin{smallmatrix}
         *               g_i + (\gamma_i/\lambda_i) \hat{e}_i \\
         *               \gamma_i/\lambda_i
         *           \end{smallmatrix}\bigr)
         */
        const value_type theta_by_lambda = -neg_theta/lambda;
        const value_type   eta_by_lambda = -neg_eta  /lambda;
        const value_type gamma_by_lambda = -neg_gamma/lambda;
        next_ehat.clear();
        next_ehat.push_back(eta_by_lambda);
        for (size_type j = 0; j < i; ++j)
        {
            s[j] += theta_by_lambda*ehat[j];
            next_ehat.push_back(ehat[j] + eta_by_lambda*g[j]);
            g[j] += gamma_by_lambda*ehat[j];
        }
        s.push_back(theta_by_lambda);
        g.push_back(gamma_by_lambda);
        ehat.swap(next_ehat);
 
        // \lambda_{i+1} = \lambda_i - \eta_i \gamma_i / \lambda_i
        lambda -= neg_eta*neg_gamma/lambda;
    }
 
    // Recursion for i = n differs slightly per Zohar's "Last Computed Values"
    // Computing g_n above was unnecessary but the incremental expense is small
    {
        reverse_iterator<RandomAccessIterator> rhat_first(r_first + n);
 
        // \theta_n =  \delta_{n+1}  - \tilde{s}_n \hat{r}_n
        const value_type neg_theta = inner_product(
                s.begin(), s.end(), rhat_first, value_type(-(*++d_first)));
 
        /*
         * s_{n+1} = \bigl(\begin{smallmatrix}
         *              s_n + (\theta_n/\lambda_n) \hat{e}_n \\
         *              \theta_n/\lambda_n
         *          \end{smallmatrix}\bigr)
         */
        const value_type theta_by_lambda = -neg_theta/lambda;
        for (size_type j = 0; j < n; ++j)
        {
            s[j] += theta_by_lambda*ehat[j];
        }
        s.push_back(theta_by_lambda);
    }
 
    // Output solution
    copy(s.begin(), s.end(), s_first);
}
 
template<class RandomAccessIterator,
         class ForwardIterator>
void zohar_linear_solve(RandomAccessIterator a_first,
                        RandomAccessIterator a_last,
                        RandomAccessIterator r_first,
                        ForwardIterator      d_first)
{
    return zohar_linear_solve(a_first, a_last, r_first, d_first, d_first);
}
 
 
template<class RandomAccessIterator,
         class ForwardIterator>
void zohar_linear_solve(RandomAccessIterator a_first,
                        RandomAccessIterator a_last,
                        ForwardIterator      d_first)
{
    return zohar_linear_solve(a_first, a_last, a_first, d_first);
}
 
 
struct mean_subtracted
{
    template <typename Result, typename Integer>
    static Result empirical_variance_zero(Integer N)
    {
        assert(N >= 1);
 
        Result den = N;
        return 1 / den;
    }
};
 
struct mean_retained
{
    template <typename Result, typename Integer>
    static Result empirical_variance_zero(Integer)
    {
        return 0;
    }
};
 
namespace { // anonymous
 
// Oh the simultaneous love/hate because -Werror => -Werror=type-limits...
 
template <typename T, bool> struct is_nonnegative_helper;
 
template <typename T> struct is_nonnegative_helper<T, /*signed?*/ true>
{
    static bool check(const T& t) { return t >= 0; }
};
 
template <typename T> struct is_nonnegative_helper<T, /*signed?*/ false>
{
    static bool check(const T&) { return true; }
};
 
template <typename T> bool is_nonnegative(const T& t)
{
    using std::numeric_limits;
    return is_nonnegative_helper<T,numeric_limits<T>::is_signed>::check(t);
}
 
}
 
 
struct estimation_method {};
 
template <class MeanHandling>
class YuleWalker : public estimation_method
{
public:
 
    template <typename Result, typename Integer1, typename Integer2>
    static Result empirical_variance(Integer1 N, Integer2 i)
    {
        assert(N >= 1);
        assert(is_nonnegative(i));
        assert(static_cast<Integer1>(i) <= N);
 
        if (i == 0)
            return MeanHandling::template empirical_variance_zero<Result>(N);
 
        Result num = N - i, den = N*(N + 2);
        return num / den;
    }
};
 
template <class MeanHandling>
class Burg : public estimation_method
{
public:
 
    template <typename Result, typename Integer1, typename Integer2>
    static Result empirical_variance(Integer1 N, Integer2 i)
    {
        assert(N >= 1);
        assert(is_nonnegative(i));
        assert(static_cast<Integer1>(i) <= N);
 
        if (i == 0)
            return MeanHandling::template empirical_variance_zero<Result>(N);
 
        Result den = N + 1 - i;
        return 1 / den;
    }
};
 
template <class MeanHandling>
class LSFB : public estimation_method
{
public:
 
    template <typename Result, typename Integer1, typename Integer2>
    static Result empirical_variance(Integer1 N, Integer2 i)
    {
        assert(N >= 1);
        assert(is_nonnegative(i));
        assert(static_cast<Integer1>(i) <= N);
 
        if (i == 0)
            return MeanHandling::template empirical_variance_zero<Result>(N);
 
        // Factorizing expression will cause problems in unsigned arithmetic
        Result den = N + Result(3)/2 - Result(3)/2 * i;
        return 1 / den;
    }
};
 
template <class MeanHandling>
class LSF : public estimation_method
{
public:
 
    template <typename Result, typename Integer1, typename Integer2>
    static Result empirical_variance(Integer1 N, Integer2 i)
    {
        assert(N >= 1);
        assert(is_nonnegative(i));
        assert(static_cast<Integer1>(i) <= N);
 
        if (i == 0)
            return MeanHandling::template empirical_variance_zero<Result>(N);
 
        // Factorizing expression will cause problems in unsigned arithmetic
        Result den = N + 2 - 2*i;
        return 1 / den;
    }
};
 
template <class EstimationMethod,
          typename Result,
          typename Integer1 = std::size_t,
          typename Integer2 = Integer1>
struct empirical_variance_function
{
    typedef Integer1 first_argument_type;
    typedef Integer2 second_argument_type;
    typedef Result   result_type;
 
    Result operator() (Integer1 N, Integer2 i) const
    {
        return EstimationMethod::template empirical_variance<Result>(N, i);
    }
};
 
template <class EstimationMethod,
          typename Result,
          typename Integer1 = std::size_t,
          typename Integer2 = Integer1>
class empirical_variance_generator
{
private:
 
    Integer1 N;
    Integer2 i;
 
public:
 
    typedef Result result_type;
 
    empirical_variance_generator(Integer1 N) : N(N), i(0) {}
 
    Result operator() ()
    {
        return EstimationMethod::template empirical_variance<Result>(N, i++);
    }
};
 
template <class EstimationMethod,
          typename Result,
          typename Integer1 = std::size_t,
          typename Integer2 = Integer1>
class empirical_variance_iterator
{
private:
    Integer1 N;
    Integer2 i;
 
public:
    typedef std::random_access_iterator_tag iterator_category;
    typedef Result                          value_type;
    typedef std::ptrdiff_t                  difference_type;
    typedef const Result*                   pointer;
    typedef const Result&                   reference;
 
    empirical_variance_iterator() : N(0), i(1) {}  // Null
 
    empirical_variance_iterator(Integer1 N) : N(N), i(0)
        { assert(N >= 1); }
 
    empirical_variance_iterator(Integer1 N, Integer2 i) : N(N), i(i) {}
 
    // Forward traversal
 
    empirical_variance_iterator& operator++()
        { assert(N >= 1); ++i; return *this; }
 
    empirical_variance_iterator operator++(int)
        { assert(N >= 1); return empirical_variance_iterator(N, i++); }
 
    empirical_variance_iterator operator+(const difference_type& k)
        { assert(N >= 1); return empirical_variance_iterator(N, i + k); }
 
    empirical_variance_iterator& operator+=(const difference_type& k)
        { assert(N >= 1); i += k; return *this; }
 
    // Backward traversal
 
    empirical_variance_iterator& operator--()
        { assert(N >= 1); --i; return *this; }
 
    empirical_variance_iterator operator--(int)
        { assert(N >= 1); return empirical_variance_iterator(N, i--); }
 
    empirical_variance_iterator operator-(const difference_type& k)
        { assert(N >= 1); return empirical_variance_iterator(N, i - k); }
 
    empirical_variance_iterator& operator-=(const difference_type& k)
        { assert(N >= 1); i -= k; return *this; }
 
    // Distance support
 
    difference_type operator-(const empirical_variance_iterator& other) const
    {
        if (!this->N) {
            return 1 + other.N - other.i;
        } else if (!other.N) {
            return -static_cast<difference_type>(1 + this->N - this->i);
        } else {
            assert(this->N == other.N);
            return this->i - other.i;
        }
    }
 
    // EqualityComparable
 
    bool operator==(const empirical_variance_iterator& other) const
    {
        if (!this->N) {
            return other.i >= static_cast<Integer2>(other.N + 1);
        } else if (!other.N) {
            return this->i >= static_cast<Integer2>(this->N + 1);
        } else {
            return this->N == other.N && this->i == other.i;
        }
    }
 
    bool operator!=(const empirical_variance_iterator& other) const
        { return !(*this == other); }
 
    // LessThanComparable will trigger assertion on nonsense N cases
 
    bool operator<(const empirical_variance_iterator& other) const
        { return (*this - other) < 0; }
 
    bool operator<=(const empirical_variance_iterator& other) const
        { return (*this - other) <= 0; }
 
    bool operator>(const empirical_variance_iterator& other) const
        { return (*this - other) > 0; }
 
    bool operator>=(const empirical_variance_iterator& other) const
        { return (*this - other) >= 0; }
 
    // Dereference operations
 
    const value_type operator*() const
    {
        assert(is_nonnegative(i));
        assert(i <= static_cast<Integer2>(N));
 
        return EstimationMethod::template empirical_variance<Result>(N, i);
    }
 
    const value_type operator[](const difference_type &k) const
    {
        assert(is_nonnegative(i + k));
        assert(i + k <= static_cast<Integer2>(N));
 
        return EstimationMethod::template empirical_variance<Result>(N, i + k);
    }
};
 
 
struct criterion
{
    template <typename Result, typename Input>
    static Result underfit_penalty(Input sigma2e)
    {
        using std::log;
        return log(Result(sigma2e));
    }
};
 
template <class    Criterion,
          typename Result,
          typename Integer1,
          typename Integer2>
Result evaluate(Result sigma2e, Integer1 N, Integer2 p)
{
    Result underfit = Criterion::template underfit_penalty<Result>(sigma2e);
    Result overfit  = Criterion::template overfit_penalty<Result>(N, p);
    return underfit + overfit;
}
 
template <int AlphaNumerator = 3, int AlphaDenominator = 1>
struct GIC : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        return Result(AlphaNumerator) * p / (N * AlphaDenominator);
    }
};
 
struct AIC : public GIC<2> {};
 
struct BIC : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        using std::log;
        return log(Result(N)) * p / N;
    }
};
 
struct MCC : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        using std::log;
        return 2*log(log(Result(N))) * p / N;
    }
};
 
struct AICC : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        return 2 * Result(p) / (N - p - 1);
    }
};
 
template <class EstimationMethod,
          int AlphaNumerator = 3,
          int AlphaDenominator = 1 >
struct FIC : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        // This accumulate invocation is inefficient but always correct
        using std::accumulate;
        typedef empirical_variance_iterator<
                EstimationMethod, Result, Integer1, Integer2
            > evi_type;
        Result sum = accumulate(evi_type(N), evi_type(N) + p, Result(0));
 
        return AlphaNumerator * sum / AlphaDenominator;
    }
};
 
template <class MeanHandling,
          int AlphaNumerator,
          int AlphaDenominator>
struct FIC<YuleWalker<MeanHandling>, AlphaNumerator, AlphaDenominator>
    : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        Result v0  = YuleWalker<MeanHandling>
            ::template empirical_variance<Result>(N, Integer2(0));
 
        Result num = (2*N - p - 1)*p;  // Avoids non-positive values
        Result den = 2*N*(N + 2);
 
        return AlphaNumerator * (v0 + num/den) / AlphaDenominator;
    }
};
 
// Specializations of the FIC for efficiency when digamma is available.
#ifdef AR_DIGAMMA
 
template <class MeanHandling,
          int AlphaNumerator,
          int AlphaDenominator>
struct FIC<Burg<MeanHandling>, AlphaNumerator, AlphaDenominator>
    : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        Result t  = Burg<MeanHandling>
            ::template empirical_variance<Result>(N, Integer2(0));
        t -= AR_DIGAMMA(N + 1);
        t += AR_DIGAMMA(N + 1 - p);
        return AlphaNumerator * t / AlphaDenominator;
    }
};
 
template <class MeanHandling,
          int AlphaNumerator,
          int AlphaDenominator>
struct FIC<LSFB<MeanHandling>, AlphaNumerator, AlphaDenominator>
    : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        Result t  = LSFB<MeanHandling>
            ::template empirical_variance<Result>(N, Integer2(0));
        Result a = AR_DIGAMMA((3 + 2*N) / Result(3)    );
        Result b = AR_DIGAMMA((3 + 2*N) / Result(3) - p);
        return AlphaNumerator * (t - 2*(a - b)/3) / AlphaDenominator;
    }
};
 
template <class MeanHandling,
          int AlphaNumerator,
          int AlphaDenominator>
struct FIC<LSF<MeanHandling>, AlphaNumerator, AlphaDenominator>
    : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        Result t  = LSF<MeanHandling>
            ::template empirical_variance<Result>(N, Integer2(0));
        Result a = AR_DIGAMMA((2 + N) / Result(2)    );
        Result b = AR_DIGAMMA((2 + N) / Result(2) - p);
        return AlphaNumerator * (t - (a - b)/2) / AlphaDenominator;
    }
};
 
#endif /* AR_DIGAMMA */
 
template <class EstimationMethod>
struct FSIC : public criterion
{
 
private:
 
    template <typename Result, typename Integer1, typename Integer2>
    class product_iterator
        : public empirical_variance_iterator<
                EstimationMethod, Result, Integer1, Integer2
          >
    {
    private:
        typedef empirical_variance_iterator<
                EstimationMethod, Result, Integer1, Integer2
            > base;
 
    public:
        typedef typename base::difference_type   difference_type;
        typedef typename base::iterator_category iterator_category;
        typedef typename base::pointer           pointer;
        typedef typename base::reference         reference;
        typedef typename base::value_type        value_type;
 
        product_iterator(Integer1 N) : base(N) {}
 
        value_type operator*() const
        {
            const value_type v = this->base::operator*();
            return (1 + v) / (1 - v);
        }
 
        value_type operator[](const difference_type &k) const
        {
            const value_type v = this->base::operator[](k);
            return (1 + v) / (1 - v);
        }
 
    };
 
public:
 
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        // This accumulate invocation is inefficient but always correct
        using std::multiplies;
        using std::accumulate;
        product_iterator<Result, Integer1, Integer2> first(N), last(N);
        last += p; // Avoids overloading requirements for (last + p)
 
        return accumulate(first, last, Result(1), multiplies<Result>()) - 1;
    }
};
 
// Specializations of the FIC for efficiency when Pochhammer is available.
#ifdef AR_POCHHAMMER
 
template <class MeanHandling>
struct FSIC<YuleWalker<MeanHandling> > : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        Result v0  = YuleWalker<MeanHandling>
            ::template empirical_variance<Result>(N, Integer2(0));
        Result a = AR_POCHHAMMER(N*(N+3) - p, p);
        Result b = AR_POCHHAMMER(Result(1 + N) - N*N, p);
 
        return (1 + v0) / (1 - v0) * (a / b) - 1;
    }
};
 
template <class MeanHandling>
struct FSIC<Burg<MeanHandling> > : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        Result v0  = Burg<MeanHandling>
            ::template empirical_variance<Result>(N, Integer2(0));
        Result a = AR_POCHHAMMER(-Result(1) - N, p);
        Result b = AR_POCHHAMMER( Result(1) - N, p);
 
        return (1 + v0) / (1 - v0) * (a / b) - 1;
    }
};
 
template <class MeanHandling>
struct FSIC<LSFB<MeanHandling> > : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        Result v0  = LSFB<MeanHandling>
            ::template empirical_variance<Result>(N, Integer2(0));
        Result a = AR_POCHHAMMER((-Result(2)-2*N)/3, p);
        Result b = AR_POCHHAMMER(( Result(2)-2*N)/3, p);
 
        return (1 + v0) / (1 - v0) * (a / b) - 1;
    }
};
 
template <class MeanHandling>
struct FSIC<LSF<MeanHandling> > : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        Result v0  = LSF<MeanHandling>
            ::template empirical_variance<Result>(N, Integer2(0));
        Result a = AR_POCHHAMMER((-Result(1)-N)/2, p);
        Result b = AR_POCHHAMMER(( Result(1)-N)/2, p);
 
        return (1 + v0) / (1 - v0) * (a / b) - 1;
    }
};
 
#endif /* AR_POCHHAMMER */
 
 
template <class EstimationMethod>
struct CIC : public criterion
{
    template <typename Result, typename Integer1, typename Integer2>
    static Result overfit_penalty(Integer1 N, Integer2 p)
    {
        using std::max;
        return max(
                FSIC<EstimationMethod>::template overfit_penalty<Result>(N, p),
                FIC <EstimationMethod>::template overfit_penalty<Result>(N, p)
            );
    }
};
 
 
template <class    Criterion,
          typename Integer1,
          typename Integer2,
          class    InputIterator,
          class    OutputIterator>
typename std::iterator_traits<InputIterator>::difference_type
evaluate_models(Integer1       N,
                Integer2       ordfirst,
                InputIterator  first,
                InputIterator  last,
                OutputIterator crit)
{
    using std::iterator_traits;
    using std::numeric_limits;
 
    typedef InputIterator iterator;
    typedef typename iterator_traits<iterator>::difference_type difference_type;
    typedef typename iterator_traits<iterator>::value_type      value_type;
 
    // Short circuit on trivial input
    if (first == last)
        return numeric_limits<difference_type>::min();
 
    // Handle first iteration without comparison as AICC blows up on N == 1
    value_type best_val = evaluate<Criterion>(*first++, N, ordfirst);
    difference_type best_pos = 0, dist = 0;
 
    // Scan through rest of candidates (up to order N-1) updating best as we go
    while (first != last && static_cast<Integer1>(++dist) < N)
    {
        value_type candidate = evaluate<Criterion>(*first++, N, dist+ordfirst);
        *crit++ = candidate;
 
        if (candidate < best_val) {
            best_val = candidate;
            best_pos = dist;
        }
    }
 
    return best_pos;
}
 
template <class    Criterion,
          typename Integer1,
          typename Integer2,
          class    Sequence1,
          class    Sequence2,
          class    Sequence3,
          class    Sequence4,
          class    OutputIterator>
typename Sequence1::difference_type
best_model(Integer1       N,
           Integer2       minorder,
           Sequence1&     params,
           Sequence2&     sigma2e,
           Sequence3&     gain,
           Sequence4&     autocor,
           OutputIterator crit)
{
    using std::advance;
    using std::copy_backward;
    using std::invalid_argument;
    using std::size_t;
 
    // Ensure all inputs have conformant sizes
    // Sanity checks written to provide order messages readable by end users
    AR_ENSURE_ARG(sigma2e.size() > 0);
    const size_t maxorder = sigma2e.size() - 1;
    AR_ENSURE_ARG(is_nonnegative(minorder));
    AR_ENSURE_ARG(static_cast<size_t>(minorder) <= maxorder);
    AR_ENSURE_ARG(params .size() == (sigma2e.size()-1)*sigma2e.size()/2);
    AR_ENSURE_ARG(gain   .size() == sigma2e.size());
    AR_ENSURE_ARG(autocor.size() == sigma2e.size());
 
    // Find the best model index according to the given minorder and criterion
    typename Sequence1::difference_type best = minorder;
    {
        typename Sequence2::iterator iter = sigma2e.begin();
        advance(iter, minorder);
        best += evaluate_models<Criterion>(N, minorder, iter,
                                           sigma2e.end(), crit);
    }
 
    // Now trim away everything and leaving only the best model in Sequences...
 
    // ...first in params...
    {
        // Best parameters might overlap beginning of params so copy_backwards.
        // AR(0) is trivial and AR(1) starts at params.begin(), hence off by 1.
        typename Sequence1::iterator first  = params.begin();
        typename Sequence1::iterator last   = params.begin();
        typename Sequence1::iterator result = params.begin();
        advance(first,  ((best-1)*best)/2       );
        advance(last,   ((best-1)*best)/2 + best);
        advance(result,                     best);
        copy_backward(first, last, result);
        params.resize(best);
    }
 
    // ...next in sigma2e...
    {
        typename Sequence2::iterator cursor = sigma2e.begin();
        advance(cursor, best);
        *sigma2e.begin() = *cursor;
        sigma2e.resize(1);
    }
 
    // ...next in gain...
    {
        typename Sequence2::iterator cursor = gain.begin();
        advance(cursor, best);
        *gain.begin() = *cursor;
        gain.resize(1);
    }
 
    // ...and last in autocor...
    autocor.resize(best + 1);
 
    // ...but notice that resizing does not free capacity for std::vectors.
 
    return best;
}
 
namespace { // anonymous
 
// Used to discard output per http://stackoverflow.com/questions/335930/
struct null_output
{
    typedef std::output_iterator_tag iterator_category;
    typedef void                     value_type;
    typedef void                     difference_type;
    typedef void                     pointer;
    typedef void                     reference;
 
    template <typename T> void operator=(const T&) {}
 
    null_output& operator++()
        { return *this; }
 
    null_output operator++(int)
        { null_output it(*this); ++*this; return it; }
 
    null_output& operator*()
        { return *this; }
};
 
}
 
template <class    Criterion,
          typename Integer1,
          typename Integer2,
          class    InputIterator>
typename std::iterator_traits<InputIterator>::difference_type
evaluate_models(Integer1      N,
                Integer2      ordfirst,
                InputIterator first,
                InputIterator last)
{
    return evaluate_models<Criterion>(N, ordfirst, first, last, null_output());
}
 
template <class    Criterion,
          typename Integer1,
          typename Integer2,
          class    Sequence1,
          class    Sequence2,
          class    Sequence3,
          class    Sequence4>
typename Sequence1::difference_type
best_model(Integer1       N,
           Integer2       minorder,
           Sequence1&     params,
           Sequence2&     sigma2e,
           Sequence3&     gain,
           Sequence4&     autocor)
{
    return best_model<Criterion>(N, minorder, params, sigma2e, gain,
                                 autocor, null_output());
}
 
template <template <class> class EstimationMethod,
          typename Integer1,
          typename Integer2,
          class    Sequence1,
          class    Sequence2 = Sequence1,
          class    Sequence3 = Sequence2,
          class    Sequence4 = Sequence3>
struct best_model_function
{
    typedef typename Sequence1::difference_type (*type)(Integer1   N,
                                                        Integer2   minorder,
                                                        Sequence1& params,
                                                        Sequence2& sigma2e,
                                                        Sequence3& gain,
                                                        Sequence4& autocor);
 
    template <class CharT, class Traits, class Allocator>
    static type lookup(std::basic_string<CharT,Traits,Allocator> abbreviation,
                       const bool subtract_mean);
 
    static type lookup(const char *abbreviation, const bool subtract_mean)
    {
        std::string s(abbreviation);
        return lookup(s, subtract_mean);
    }
};
 
template <template <class> class EstimationMethod,
          typename Integer1,
          typename Integer2,
          class    Sequence1,
          class    Sequence2,
          class    Sequence3,
          class    Sequence4>
template <class CharT, class Traits, class Allocator>
typename best_model_function<
    EstimationMethod,Integer1,Integer2,Sequence1,Sequence2,Sequence3,Sequence4
>::type
best_model_function<
    EstimationMethod,Integer1,Integer2,Sequence1,Sequence2,Sequence3,Sequence4
>::lookup(std::basic_string<CharT,Traits,Allocator> abbrev,
          const bool subtract_mean)
{
    using std::basic_string;
    using std::toupper;
 
    // Canonicalize the abbreviation by making it uppercase and trimming it
    // For nothing but this reason the method accepts 'abbrev' by value
    typedef basic_string<CharT,Traits,Allocator> string_type;
    for (typename string_type::iterator p = abbrev.begin();
         abbrev.end() != p; ++p)
    {
        *p = toupper(*p);
    }
    abbrev.erase(0, abbrev.find_first_not_of(" \n\r\t"));
    abbrev.erase(1 + abbrev.find_last_not_of(" \n\r\t"));
 
    // Obtain best_model(...) per abbrev, subtract_mean, EstimationMethod
    // and assign to retval to provide unambiguous overload resolution.
    type retval;
 
    if      (abbrev.empty() || 0 == abbrev.compare("CIC" ))  // Default
    {
        if (subtract_mean)
            retval = best_model<CIC<EstimationMethod<mean_subtracted> > >;
        else
            retval = best_model<CIC<EstimationMethod<mean_retained  > > >;
    }
    else if (0 == abbrev.compare("AIC" ))
    {
        retval = best_model<AIC>;
    }
    else if (0 == abbrev.compare("AICC"))
    {
        retval = best_model<AICC>;
    }
    else if (0 == abbrev.compare("BIC" ))
    {
        retval = best_model<BIC>;
    }
    else if (0 == abbrev.compare("FIC" ))
    {
        if (subtract_mean)
            retval = best_model<FIC<EstimationMethod<mean_subtracted> > >;
        else
            retval = best_model<FIC<EstimationMethod<mean_retained  > > >;
    }
    else if (0 == abbrev.compare("FSIC"))
    {
        if (subtract_mean)
            retval = best_model<FSIC<EstimationMethod<mean_subtracted> > >;
        else
            retval = best_model<FSIC<EstimationMethod<mean_retained  > > >;
    }
    else if (0 == abbrev.compare("GIC" ))
    {
        retval = best_model<GIC<> >;
    }
    else if (0 == abbrev.compare("MCC" ))
    {
        retval = best_model<MCC>;
    } else
    {
        retval = NULL;
    }
 
    return retval;
}
 
template<class Iterator>
class strided_adaptor
{
public:
    typedef typename std::iterator_traits<Iterator>::iterator_category iterator_category;
    typedef typename std::iterator_traits<Iterator>::value_type        value_type;
    typedef typename std::iterator_traits<Iterator>::difference_type   difference_type;
    typedef typename std::iterator_traits<Iterator>::pointer           pointer;
    typedef typename std::iterator_traits<Iterator>::reference         reference;
 
    strided_adaptor() : i(NULL), step(0) {};
 
    strided_adaptor(const strided_adaptor& o) : i(o.i), step(o.step) {}
 
    strided_adaptor(Iterator x, difference_type n) : i(x), step(n) {}
 
    strided_adaptor& operator++()
    {
        using std::advance;
        advance(i, step);
        return *this;
    }
 
    strided_adaptor operator++(int)
    {
        strided_adaptor tmp(*this);
        using std::advance;
        advance(i, step);
        return tmp;
    }
 
    strided_adaptor& operator+=(difference_type x)
    {
        using std::advance;
        advance(i, x*step);
        return *this;
    }
 
    strided_adaptor& operator--( )
    {
        using std::advance;
        advance(i, -step);
        return *this;
    }
 
    strided_adaptor operator--(int)
    {
        strided_adaptor tmp(*this);
        using std::advance;
        advance(i, -step);
        return tmp;
    }
 
    strided_adaptor& operator-=(difference_type x)
    {
        using std::advance;
        advance(i, -x*step);
        return *this;
    }
 
    reference operator[](difference_type n)
    {
        return i[n * step];
    }
 
    reference operator*()
    {
        return *i;
    }
 
    bool operator==(const strided_adaptor& o) const
    {
        return i == o.i;
    }
 
    bool operator!=(const strided_adaptor& o) const
    {
        return i != o.i;
    }
 
    bool operator<(const strided_adaptor& o) const
    {
        return i < o.i;
    }
 
    difference_type operator-(const strided_adaptor& o) const
    {
        return (i - o.i) / step;
    }
 
    strided_adaptor operator+(difference_type x) const
    {
        strided_adaptor tmp(*this);
        tmp += x;
        return tmp;
    }
 
private:
    Iterator i;
    difference_type step;
};
 
 
} // end namespace ar
 
#endif /* AR_HPP */