Multithreaded Monte-Carlo Integration

Question

I have a PR implementing multithreaded naive Monte-Carlo integration here. My requirements for the class are the following:

1. It should support progress reporting, ETA, and graceful cancellation.
2. It should handle finite and infinite intervals.
3. It should allow restarts, and changing of error goal in flight.

Unfortunately, the imposition of these requirements has left me dissatisfied with the design and look of the code:

#ifndef BOOST_MATH_QUADRATURE_NAIVE_MONTE_CARLO_HPP
#define BOOST_MATH_QUADRATURE_NAIVE_MONTE_CARLO_HPP
#include <algorithm>
#include <vector>
#include <atomic>
#include <functional>
#include <future>
#include <thread>
#include <initializer_list>
#include <utility>
#include <random>
#include <chrono>
#include <map>

namespace boost { namespace math { namespace quadrature {

namespace detail {
  enum class limit_classification {FINITE,
                                   LOWER_BOUND_INFINITE,
                                   UPPER_BOUND_INFINITE,
                                   DOUBLE_INFINITE};
}

template<class Real, class F, class Policy = boost::math::policies::policy<>>
class naive_monte_carlo
{
public:
    naive_monte_carlo(const F& f,
                      std::vector<std::pair<Real, Real>> const & bounds,
                      Real error_goal,
                      size_t threads = std::thread::hardware_concurrency()): m_num_threads{threads}
    {
        using std::isfinite;
        using std::numeric_limits;
        size_t n = bounds.size();
        m_lbs.resize(n);
        m_dxs.resize(n);
        m_limit_types.resize(n);
        m_volume = 1;
        static const char* function = "boost::math::quadrature::naive_monte_carlo<%1%>";
        for (size_t i = 0; i < n; ++i)
        {
            if (bounds[i].second <= bounds[i].first)
            {
                boost::math::policies::raise_domain_error(function, "The upper bound is <= the lower bound.\n", bounds[i].second, Policy());
                return;
            }
            if (bounds[i].first == -numeric_limits<Real>::infinity())
            {
                if (bounds[i].second == numeric_limits<Real>::infinity())
                {
                    m_limit_types[i] = detail::limit_classification::DOUBLE_INFINITE;
                }
                else
                {
                    m_limit_types[i] = detail::limit_classification::LOWER_BOUND_INFINITE;
                    // Ok ok this is bad:
                    m_lbs[i] = bounds[i].second;
                    m_dxs[i] = std::numeric_limits<Real>::quiet_NaN();
                }
            }
            else if (bounds[i].second == numeric_limits<Real>::infinity())
            {
                m_limit_types[i] = detail::limit_classification::UPPER_BOUND_INFINITE;
                m_lbs[i] = bounds[i].first;
                m_dxs[i] = std::numeric_limits<Real>::quiet_NaN();
            }
            else
            {
                m_limit_types[i] = detail::limit_classification::FINITE;
                m_lbs[i] = bounds[i].first;
                m_dxs[i] = bounds[i].second - m_lbs[i];
                m_volume *= m_dxs[i];
            }
        }

        m_f = [this, &f](std::vector<Real> & x)->Real
        {
            Real coeff = m_volume;
            for (size_t i = 0; i < x.size(); ++i)
            {
                // Variable transformation are listed at:
                // https://en.wikipedia.org/wiki/Numerical_integration
                if (m_limit_types[i] == detail::limit_classification::FINITE)
                {
                    x[i] = m_lbs[i] + x[i]*m_dxs[i];
                }
                else if (m_limit_types[i] == detail::limit_classification::UPPER_BOUND_INFINITE)
                {
                    Real t = x[i];
                    Real z = 1/(1-t);
                    coeff *= (z*z);
                    x[i] = m_lbs[i] + t*z;
                }
                else if (m_limit_types[i] == detail::limit_classification::LOWER_BOUND_INFINITE)
                {
                    Real t = x[i];
                    Real z = 1/t;
                    coeff *= (z*z);
                    x[i] = m_lbs[i] + (t-1)*z;
                }
                else
                {
                    Real t = 2*x[i] - 1;
                    Real tsq = t*t;
                    Real z = 1/(1-t);
                    z /= (1+t);
                    coeff *= 2*(1+tsq)*z*z;
                    x[i] = t*z;
                }
            }
            return coeff*f(x);
        };

        // If we don't do a single function call in the constructor,
        // we can't do a restart.
        std::vector<Real> x(m_lbs.size());
        std::random_device rd;
        std::mt19937_64 gen(rd());
        Real inv_denom = (Real) 1/( (Real) gen.max() + (Real) 2);

        if (m_num_threads == 0)
        {
            m_num_threads = 1;
        }
        Real avg = 0;
        for (size_t i = 0; i < m_num_threads; ++i)
        {
            for (size_t j = 0; j < m_lbs.size(); ++j)
            {
                x[j] = (gen()+1)*inv_denom;
                while (x[j] < std::numeric_limits<Real>::epsilon() || x[j] > 1 - std::numeric_limits<Real>::epsilon())
                {
                    x[j] = (gen()+1)*inv_denom;
                }
            }
            Real y = m_f(x);
            m_thread_averages.emplace(i, y);
            m_thread_calls.emplace(i, 1);
            m_thread_Ss.emplace(i, 0);
            avg += y;
        }
        avg /= m_num_threads;
        m_avg = avg;

        m_error_goal = error_goal;
        m_start = std::chrono::system_clock::now();
        m_done = false;
        m_total_calls = m_num_threads;
        m_variance = std::numeric_limits<Real>::max();
    }

    std::future<Real> integrate()
    {
        // Set done to false in case we wish to restart:
        m_done = false;
        return std::async(std::launch::async,
                          &naive_monte_carlo::m_integrate, this);
    }

    void cancel()
    {
        m_done = true;
    }

    Real variance() const
    {
        return m_variance.load();
    }

    Real current_error_estimate() const
    {
        using std::sqrt;
        return sqrt(m_variance.load()/m_total_calls.load());
    }

    std::chrono::duration<Real> estimated_time_to_completion() const
    {
        auto now = std::chrono::system_clock::now();
        std::chrono::duration<Real> elapsed_seconds = now - m_start;
        Real r = this->current_error_estimate()/m_error_goal.load();
        if (r*r <= 1) {
            return 0*elapsed_seconds;
        }
        return (r*r - 1)*elapsed_seconds;
    }

    void update_target_error(Real new_target_error)
    {
        m_error_goal = new_target_error;
    }

    Real progress() const
    {
        Real r = m_error_goal.load()/this->current_error_estimate();
        if (r*r >= 1)
        {
            return 1;
        }
        return r*r;
    }

    Real current_estimate() const
    {
        return m_avg.load();
    }

    size_t calls() const
    {
        return m_total_calls.load();
    }

private:

    Real m_integrate()
    {
        m_start = std::chrono::system_clock::now();
        std::vector<std::thread> threads(m_num_threads);
        for (size_t i = 0; i < threads.size(); ++i)
        {
            threads[i] = std::thread(&naive_monte_carlo::m_thread_monte, this, i);
        }
        do {
            std::this_thread::sleep_for(std::chrono::milliseconds(500));
            size_t total_calls = 0;
            for (size_t i = 0; i < m_num_threads; ++i)
            {
                total_calls += m_thread_calls[i];
            }
            Real variance = 0;
            Real avg = 0;
            for (size_t i = 0; i < m_num_threads; ++i)
            {
                size_t t_calls = m_thread_calls[i];
                // Will this overflow? Not hard to remove . . .
                avg += m_thread_averages[i]*( (Real) t_calls/ (Real) total_calls);
                variance += m_thread_Ss[i];
            }
            m_avg = avg;
            m_variance = variance/(total_calls - 1);
            m_total_calls = total_calls;
            // Allow cancellation:
            if (m_done)
            {
                break;
            }
        } while (this->current_error_estimate() > m_error_goal);
        // Error bound met; signal the threads:
        m_done = true;
        std::for_each(threads.begin(), threads.end(),
                      std::mem_fn(&std::thread::join));
        if (m_exception)
        {
            std::rethrow_exception(m_exception);
        }
        // Incorporate their work into the final estimate:
        size_t total_calls = 0;
        for (size_t i = 0; i < m_num_threads; ++i)
        {
            total_calls += m_thread_calls[i];
        }
        Real variance = 0;
        Real avg = 0;
        for (size_t i = 0; i < m_num_threads; ++i)
        {
            size_t t_calls = m_thread_calls[i];
            // Will this overflow? Not hard to remove . . .
            avg += m_thread_averages[i]*( (Real) t_calls/ (Real) total_calls);
            variance += m_thread_Ss[i];
        }
        m_avg = avg;
        m_variance = variance/(total_calls - 1);
        m_total_calls = total_calls;

        return m_avg.load();
    }

    void m_thread_monte(size_t thread_index)
    {
        try
        {
            std::vector<Real> x(m_lbs.size());
            std::random_device rd;
            // Should we do something different if we have no entropy?
            // Apple LLVM version 9.0.0 (clang-900.0.38) has no entropy,
            // but rd() returns a reasonable random sequence.
            // if (rd.entropy() == 0)
            // {
            //     std::cout << "OMG! we have no entropy.\n";
            // }
            std::mt19937_64 gen(rd());
            Real inv_denom = (Real) 1/( (Real) gen.max() + (Real) 2);
            Real M1 = m_thread_averages[thread_index];
            Real S = m_thread_Ss[thread_index];
            // Kahan summation is required. See the implementation discussion.
            Real compensator = 0;
            size_t k = m_thread_calls[thread_index];
            while (!m_done)
            {
                int j = 0;
                // If we don't have a certain number of calls before an update, we can easily terminate prematurely
                // because the variance estimate is way too low.
                int magic_calls_before_update = 2048;
                while (j++ < magic_calls_before_update)
                {
                    for (size_t i = 0; i < m_lbs.size(); ++i)
                    {
                        x[i] = (gen()+1)*inv_denom;
                        while (x[i] < std::numeric_limits<Real>::epsilon() || x[i] > 1 - std::numeric_limits<Real>::epsilon())
                        {
                            x[i] = (gen()+1)*inv_denom;
                        }
                    }
                    Real f = m_f(x);
                    ++k;
                    Real term = (f - M1)/k;
                    Real y1 = term - compensator;
                    Real M2 = M1 + y1;
                    compensator = (M2 - M1) - y1;
                    S += (f - M1)*(f - M2);
                    M1 = M2;
                }
                m_thread_averages[thread_index] = M1;
                m_thread_Ss[thread_index] = S;
                m_thread_calls[thread_index] = k;
            }
        }
        catch (...)
        {
            // Signal the other threads that the computation is ruined:
            m_done = true;
            m_exception = std::current_exception();
        }
    }

    std::function<Real(std::vector<Real> &)> m_f;
    size_t m_num_threads;
    std::atomic<Real> m_error_goal;
    std::atomic<bool> m_done;
    std::vector<Real> m_lbs;
    std::vector<Real> m_dxs;
    std::vector<detail::limit_classification> m_limit_types;
    Real m_volume;
    std::atomic<size_t> m_total_calls;
    // I wanted these to be vectors rather than maps,
    // but you can't resize a vector of atomics.
    std::map<size_t, std::atomic<size_t>> m_thread_calls;
    std::atomic<Real> m_variance;
    std::map<size_t, std::atomic<Real>> m_thread_Ss;
    std::atomic<Real> m_avg;
    std::map<size_t, std::atomic<Real>> m_thread_averages;
    std::chrono::time_point<std::chrono::system_clock> m_start;
    std::exception_ptr m_exception;
};

}}}
#endif

I have the following specific complaints:

1. I feel the classification of the limits as finite, half-infinite, and infinite is an unnatural hack. Is this necessary?
2. Each thread needs to accumulate a variance, and average, and a number of calls. These must be atomic so they can be reduced by the master thread without a race condition. However, a vector of atomics cannot be resized, so I used a map, which, though not a catastrophe, seems suboptimal. Is there a workaround?
3. I'm using (say) std::atomic<double>, which seems to have widespread compiler support, but won't have official status until C++20. A workaround is to use a mutex, but a mutex is a disaster for performance. What should be done?
4. I'm taking the function by const &, but should it be forwarded &&? Or should it provide two constructors?

5. I tried many random number generators, and the Mersenne twister seems to be the only one that doesn't contract 'weird seed' and reliably 'does the right thing'. However, it returns a 64-bit integer, which must be remapped into the open interval ]0, 1[ (or else it'll hit singularities on the boundary). Using std::uniform_real_distribution was too slow, so I use (gen()+gen.min()+1)/(gen.max()+gen.min() + 2), which is always <1 and >0 in double precision. In float precision, it gets rounded to 1 or 0 quite frequently, so I've added the following hack:

   while (x[i] < numeric_limits<Real>::epsilon() || x[i] > 1 - numeric_limits<Real>::epsilon()) { x[i] = (gen()+1)*inv_denom;}

   Yuck.

6. Generally, it is bad form to take STL containers as arguments, as the bounds require. Any suggestions for replacement?

7. Perhaps the template parameter Real is redundant, and can be replaced by the return type of the function?

Answer 1

Since C++17 (and its now 2018) You don;t need to do this ugly thing

namespace boost { namespace math { namespace quadrature {
}}}

You can simply do this:

namespace boost::math::quadrature
{
    int x = 5;
}

Personal thing I don't like this:

    using std::isfinite;
    using std::numeric_limits;

But you did it perfectly by enclosing it inside the smallest scope possible (so I can really say anything bad about it and its absolutely fine in this context). My personal preference is that anything from the namespace/class that I am currently defining does not have a prefix. But Anything from another namespace needs a prefix (though a namespace alias is fine for long things I like the prefix to be short 3 or 4 letters).

    namespace asio = boost::asio;  // example

Note:

This is never used:

    using std::isfinite;

This is used inconsistently:

    using std::numeric_limits;


        if (bounds[i].first == -numeric_limits<Real>::infinity())
                             // ^^^ No std:: here
        m_dxs[i] = std::numeric_limits<Real>::quiet_NaN();
          //       ^^^^^ but here you use it.

These two could have been done as part of the initializer list

    // Doing this is potentially slightly less efficient.
    // But compiler will probably spot that. So not a huge thing.
    m_lbs.resize(n);
    m_dxs.resize(n);
    m_limit_types.resize(n);

But I like using the initializer list whenever possible.

Don't like this:

    static const char* function = "boost::math::quadrature::naive_monte_carlo<%1%>";

Especially since it is only used once. Also most compilers have a macro that does something similar automatically.

Really don't like this:

   m_f = [this, &f](std::vector<Real> & x)->Real

That just makes things hard to read. Store f as the member and change m_f into a normal method that uses f (but have a better name). I was looking through to find where f was used you know how many false positives I need to check before I found the point were you called the function (put a meaningful name on the function userFunc?).

You use a random device in several places:

    std::vector<Real> x(m_lbs.size());
    std::random_device rd;
    std::mt19937_64 gen(rd());

Each time you create a new one. Best to only create this once (it is relatively expensive to create). I would even inject the random device into the object rather than let the device create its own (then you can control the randomness for testing).

Easier to read this

    if (m_num_threads == 0)
    {
        m_num_threads = 1;
    }

When it is written as:

    m_num_threads = std::max(m_num_threads, 1);

Answer 2

A few observations.

One possibility to avoid the limits classification (concern #1) would be to use 4 different vectors to hold the values for the 4 classifications. Then, rather than having a series of if statements in your m_f loop, you'd have 4 loops (one for each limit classification).

Your loops to generate your random numbers (concern #5) can avoid the duplicate code by using a do/while loop instead of a while:

do {
    x[j] = (gen()+1)*inv_denom;
} while (x[j] < std::numeric_limits<Real>::epsilon() || x[j] > 1 - std::numeric_limits<Real>::epsilon());

This could also be pushed off into a function to be called in both places (pass in a reference to x and gen).

Each thread should have its own random number generator.

Your usage of maps for the thread variables (concern #2) is valid. In its current state, you still have a race condition, because a thread can be in the midst of updating those three variables (m_thread_averages, m_thread_Ss, and m_thread_calls) while the main thread in m_integrate is reading them (in a different order, m_thread_calls, m_thread_averages, then m_thread_Ss) so there is a slight possibility of getting updated values for some of those, and stale ones for others.

In m_thread_monte, each thread (instance) can keep a pointer or reference to the three m_thread variables it updates, so it won't have to look them up in the map every time. This will reduce the time taken to update these three values, reducing the race possibility mentioned above.

Answer 3

You should not write the Kahan summation algorithm inline. Rather define it as a template class and use it like this:

auto sum = kahan<double>();
sum += 1.0;
sum += 1.0e-30;
sum += 31415.0;
std::cout << sum << "\n";

This will make the rest of the code easier to read and raises it to a level of abstraction where you describe the mathematical properties instead of implementation details.