Boost math (ibeta_inv function) not thread safe?

Question

I have compiled part of boost - the ibeta_inv function - into a .Net 64 bit assembly and it worked great until I started calling it from multiple threads. Then it occationally return wrong results.

I complied it using this code (C++/CLI):

// Boost.h

#pragma once

#include <boost/math/special_functions/beta.hpp>

using namespace boost::math;

namespace Boost {

    public ref class BoostMath
    {
    public:
        double static InverseIncompleteBeta( double a, double b, double x )
        {
            return ibeta_inv(a,b,x);
        }
    };
}

Has anyone tried this before?

I have NOT tried this outside .Net, so I don't know if this is the cause, but I really don't see why, since it works great single threaded.

Usage (C#):

private void calcBoost(List<Val> vals)
{
    //gives WRONG results (sometimes):
    vals.AsParallel().ForAll(v => v.BoostResult = BoostMath.InverseIncompleteBeta(v.A, v.B, v.X));
    //gives CORRECT results:
    vals.ForEach(v => v.BoostResult = BoostMath.InverseIncompleteBeta(v.A, v.B, v.X));
}

UPDATE: As can be seen in my comments below - I'm not sure at all anymore that this is a Boost problem. Maybe it is some weird PLinq to C++/CLI bug??? I'm blaffed and will return with more facts later.

Growler

Documentation says whole boost.math should be thread safe as long as you use built-in floating point types (which, as I see, you do). Maybe you should file a bug? boost.org/doc/libs/release/libs/math/doc/sf_and_dist/html/… – Quackenbush

If nothing else turns up, I may try it out in a native C++ app and see if the problem remains. If so, a bug report may be the only thing to do, since I can't spot anything in the source code. Pity though... it runs twice as fast as our current implementation of the inverse incomplete beta function. – Apocynaceous

Interesting! Two thoughts come to mind: (1) treble-check you've built boost in multithreaded mode (current versions still make the distinction), and (2) this quote, from the documentation @Quackenbush linked: The reason for the latter limitation is the need to initialise symbolic constants using constructs ... as in this case there is a need for T's constructors to be run, leading to potential race conditions. I openly wonder if your code exposes this race condition unexpectedly, and I wholeheartedly back stanwise on reporting this as a bug. – Buonaparte

Lets see your Val class. Or try it at least with an immutable Val class. – Mattias

@weston: you're on to something! my Val class just contains 5 doubles. Nothing else. However, your inquiry caused me to change from class to stuct: and behold - I can't reproduce the problem when Val is a struct!!! This, however raises a bunch of new questions. I hope I can find time within a day to update this question to reflect my new findings. – Apocynaceous

Another thought: If it's not the boost library, you should be able to replace the call to the library with some simple equation and still see the problem I would assume. That could rule out boost at least. – Mattias

Exactly - that's what I've been doing - comparing to a library called alglib and I never seen any problems there... But right now all of my testing results are so weird that I'm starting to think that there must be more than one problem lurking behind the scenes... (things going wrong approx 1 time out of 10 millions rounds...) – Apocynaceous

Hi! Is there any kind of regularity in wrong results? Like every second one is wrong? It might help us to think in correct direction if you will be able to trace out which thread is returning wrong results. I had an issue with PLINQ and found out that from all 4 threads only main thread was behaving as expected. It helped me to find the bug. See my PLINQ problem here You can very number of concurrent threads by calling `.AsParallel().WithDegreeOfParalelism(n)'. Try also n=1. – Land

It is possible that this question is based on false assumptions and that the definition of "correct", which comes from well tested production code, may still not be correct!!! When test results keeps getting weirder, one must take a step back and question what the heck is going on. :-( – Apocynaceous

If you are at a loss still, you could test your RAM. A long time back I had a weird problem with CVS where a character was getting changed occasionally. A guy suggested we ran a RAM test, I'm sorry I can't remember which one he used, but there's a few out there and sure enough it showed a fault. We changed RAM and everything was good again. It's a long shot but it does happen, you could also try your test on another machine to rule out RAM and more. Windows 7 has one built in: windows.microsoft.com/en-US/windows7/… – Mattias

Answer 1

I happen to have encapsulated part of boost in a C++/CLI 64bit project and use it from C# exactly as you do.

So I threw in your C++ class in my own Boost wrapper and added this code to the C# project:

    private class Val
    {
        public double A;
        public double B;
        public double X;
        public double ParallellResult;
    }

    private static void findParallelError()
    {
        var r = new Random();

        while (true)
        {
            var vals = new List<Val>();
            for (var i = 0; i < 1000*1000; i++)
            {
                var val = new Val();
                val.A = r.NextDouble()*100;
                val.B = val.A + r.NextDouble()*1000;
                val.X = r.NextDouble();
                vals.Add(val);
            }

            // parallel calculation
            vals.AsParallel().ForAll(v => v.ParallellResult = Boost.BoostMath.InverseIncompleteBeta(v.A, v.B, v.X));

            /sequential verification
            var error = vals.Exists(v => v.ParallellResult != Boost.BoostMath.InverseIncompleteBeta(v.A, v.B, v.X));
            if (error)
                return;
        }
    }

And it just executes "forever". The parallel results are equal to the sequential results all the time. No thread unsafety here...

May I suggest that you download a fresh copy of Boost and include it in a completely new project and try this out?

I also notice that you called your result "BoostResult"... and mention in the comments something about "our current implementaion". Exactly what are you comparing your result agains? What is your definition of "correct"?

Eckart

Bonkers. I'll have to write up an apology for this question. It is based on completely wrong assumptions. A tribute to the powerfulness of Murphys Law one could say. The error is in the production wcode which I was trying to speed up. – Apocynaceous

Answer 2

It is vital that the Val class is threadsafe.

A simple way to ensure this would be to make it immutable, but I see you also have BoostResult that needs to be written. So that will need to be volatile, or have some form of locking.

public sealed class Val
{
    // Immutable fields are inheriently threadsafe 
    public readonly double A;
    public readonly double B;
    public readonly double X;

    // volatile is an easy way to make a single field thread safe
    // box and unbox to allow double as volatile
    private volatile object boostResult = 0.0;

    public Val(double A, double B, double X)
    {
        this.A = A;
        this.B = B;
        this.X = X;
    }

    public double BoostResult
    {
        get
        {
            return (double)boostResult;
        }
        set
        {
            boostResult = value;
        }
    }
}

Lock version: (see this question to determine which is best suited for your application)

public sealed class Val
{
    public readonly double A;
    public readonly double B;
    public readonly double X;

    private readonly object lockObject = new object();
    private double boostResult;

    public Val(double A, double B, double X)
    {
        this.A = A;
        this.B = B;
        this.X = X;
    }

    public double BoostResult
    {
        get
        {
            lock (lockObject)
            {
                return boostResult;
            }
        }
        set
        {
            lock (lockObject)
            {
                boostResult = value;
            }
        }
    }
}

And if you think 6 million locks will be slow, just try this:

using System;

namespace ConsoleApplication17
{
    class Program
    {
        static void Main(string[] args)
        {
            { //without locks
                var startTime = DateTime.Now;
                int i2=0;
                for (int i = 0; i < 6000000; i++)
                {
                    i2++;
                }
                Console.WriteLine(i2);
                Console.WriteLine(DateTime.Now.Subtract(startTime)); //0.01 seconds on my machine
            }
            { //with locks
                var startTime = DateTime.Now;
                var obj = new Object();
                int i2=0;
                for (int i = 0; i < 6000000; i++)
                {
                    lock (obj)
                    {
                        i2++;
                    }
                }
                Console.WriteLine(i2);
                Console.WriteLine(DateTime.Now.Subtract(startTime)); //0.14 seconds on my machine, and this isn't even in parallel.
            }
            Console.ReadLine();
        }
    }
}

Mattias

No, this is not correct. With version one - I still get the same problem (and a double cannot be volatile, but that is irrelevant since I don't read from it until all calculations are done). Version two is out of the question since I'm making this call six million times and if I'm going to lock each time - our current implementation will be much faster. BUT: making Val a struct WORKS!!! – Apocynaceous

@danbystrom Oh yeah, sorry. Have updated volatile version. Unless you are happy with your struct version, I would be interested if the new volatile version works. You may not read it, but if it's written by the main thread, even just during the initalization of a new Val, it could be cached by that thread I believe. – Mattias

I don't think you got this volatile concept correct. It tells the compiler that the compiler can't produce code that caches the variable in a register - since it can change any time. It has no more profound "thread safe" meaing. Not to my knowledge at least... Feel free to correct me! And, no, I'm not happy with my struct solution... until I know what the bonkers is going on!!! But it has encourage me to believe that I can get this working... somehow! Not sure I would have tried the struct thing if it weren't for your comment! – Apocynaceous

@danbystrom Yes, certainly caching in a register is one optimisation volatile stops. But each .net thread should be thought of to have it's own cache of memory. This is a good article. – Mattias

@danbystrom I just tried 6 million locks, in series, not even in parallel. It tool just 0.14 seconds. So I don't imagine that will dent your performance. Certainly not to the extent that parallel has no benefit. – Mattias