8
\$\begingroup\$

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?

\$\endgroup\$
  • \$\begingroup\$ Possible borderline bug: Is there any circumstance under which: return 0*elapsed_seconds; can return a non-zero value? \$\endgroup\$ – Jerry Coffin Jan 12 '18 at 19:27
  • \$\begingroup\$ @JerryCoffin: My unit tests didn't pick that up, but you're right that 0*elapsed_seconds is ugly. I did that to preserve the return type using the algebra defined on std::chrono time points. \$\endgroup\$ – user14717 Jan 12 '18 at 19:33
  • \$\begingroup\$ Hmm...okay. I'll have to think a bit more about whether there's a cleaner way under the circumstances. \$\endgroup\$ – Jerry Coffin Jan 12 '18 at 19:49
5
\$\begingroup\$

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);
\$\endgroup\$
  • \$\begingroup\$ Changed f to integrand, which is way more precise. Great suggestion. Using a single RNG is a good idea, but I saw no guarantees that the operator()() was threadsafe in the docs. Not sure that the cleaner namespace syntax is worth the trade of requiring a C++17 capable compiler, but we'll see what the maintainers say. \$\endgroup\$ – user14717 Jan 12 '18 at 22:29
1
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Used the do-while, but ideally I'd still like to only generate random numbers in the acceptable interval, rather than generate numbers and test that they are in the acceptable range. As to the race condition, since it's atomic it never sees an invalid value, though the value it sees does depend on timing. This is benign, and allows for one thread to be significantly slower than the other. \$\endgroup\$ – user14717 Jan 13 '18 at 4:58
  • 1
    \$\begingroup\$ @user14717 I was thinking you could try gen() * (1 - 2 * std::numeric_limits<Real>::epsilon()) + std::numeric_limits<Real>::epsilon() to generate the random value, but since that may be calculated using double precision (even for floats) the rounding issue might still exist. This would also skew the distribution a bit, since some numbers would be more probable than others. Since the out-of-range generated value is relatively rare, using the loop and regenerating should not be much of an issue. \$\endgroup\$ – 1201ProgramAlarm Jan 13 '18 at 6:05
  • \$\begingroup\$ @user14717 The race condition comes about because m_thread_averages can have a value from 2 loop iterations (4096 generated numbers), and m_thread_calls could have the value from only 1 loop (2048 generated numbers), which can cause the computed avg and variance to be incorrect. \$\endgroup\$ – 1201ProgramAlarm Jan 13 '18 at 6:07
  • \$\begingroup\$ Ah yes I see. But that only affects the progress reporting, not the final value. Still might be worth trying to clean up. \$\endgroup\$ – user14717 Jan 14 '18 at 15:32
1
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Great suggestion. I worry, however, that if I put a Kahan class into boost there will be an expectation that operator+= will be threadsafe. In that case, there's no point in each individual thread computing it's own average . . . well, yeah that might be way better, but there might be a lot of contention on the atomic variable. Need to perf test that one. \$\endgroup\$ – user14717 Jan 13 '18 at 15:21
  • \$\begingroup\$ Since nobody expects std::vector to be threadsafe, why should the kahan class be? \$\endgroup\$ – Roland Illig Jan 14 '18 at 7:57
  • \$\begingroup\$ You could also mark the kahan class internal (maybe by putting it in a namespace called internal); there are surely people who did something like this before. \$\endgroup\$ – Roland Illig Jan 14 '18 at 8:06
  • \$\begingroup\$ You're right: boost.org/doc/libs/1_50_0/boost/accumulators/statistics/… \$\endgroup\$ – user14717 Jan 14 '18 at 16:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.