Deterministic multithreading usage of itkMersenneTwisterRandomVariate…

[michalmeszaroswork]

…Generator.cxx

When using ITK library along with SuperElastix
https://github.com/SuperElastix/elastix
by running multiple elastix::ELASTIX objects in different threads such as :

const auto runElastix = {
elastix::ELASTIX elastixFilter;
int error = elastixFilter.RegisterImages(
fixedImg,
movingImg,
propertiesMap,
outDirPath,
false,
false,
nullptr,
nullptr);
}

std::vector<std::thread> threads;
for (int i = 0; i <= 1; ++i) {
    threads.emplace_back(runElastix);
}

for (auto &t : threads) {
    t.join();
}

they share the same instance of MersenneTwisterRandomVariateGenerator, which cause that results are nondeterministic. This simple fix stabilize results along multiple threads and does not seems to break anything in ITK.

PR Checklist

• No API changes were made (or the changes have been approved)
• No major design changes were made (or the changes have been approved)
• Added test (or behavior not changed)
• Updated API documentation (or API not changed)
• Added license to new files (if any)
• Added Python wrapping to new files (if any) as described in ITK Software Guide Section 9.5
• Added ITK examples for all new major features (if any)

Refer to the ITK Software Guide for
further development details if necessary.

[github-actions]

Thank you for contributing a pull request! 🙏

Welcome to the ITK community! 🤗👋☀️

We are glad you are here and appreciate your contribution. Please keep in mind our community participation guidelines. 📜
More support and guidance on the contribution process can be found in our contributing guide. 📖

This is an automatic message. Allow for time for the ITK community to be able to read the pull request and comment
on it.

[dzenanz]

This breaks the singleton mechanics, e.g. setting a seed from one thread will now only apply that seed to that thread. Niels, does this change break any tests in Elastix?

[N-Dekker]

Interesting, thanks @michalmeszaroswork! I just gave it a try: https://github.com/SuperElastix/elastix/tree/michalmeszaroswork-patch-1 (should take an hour or so to pass the test).

[blowekamp]

Nice job identifying the dependency onto global random. However, change the behavior of this may not be desirable by all even if it does not immediately "break" things. Let us consider another option.

The alternative is that each algorithm initializes its own random number generate ( either seeded by a user provided value, clock or the next seed from the global ). Using an algorithm specific generator removes the dependency on the global ( even a thread local is still a global and could cause issues is some cases, especially with certain thread pool or task based threading models.)

This breaks the singleton mechanics, e.g. setting a seed from one thread will now only apply that seed to that thread. Niels, does this change break any tests in Elastix?

[N-Dekker]

Nitpick: I would name it threadLocalInstance, rather than threadSafeInstance. The original instance m_StaticInstance was thread-safe already, but your essential change is to make it thread local.

Then I would suggest to have it initialized during its declaration, so that it can be const:

  const thread_local MersenneTwisterRandomVariateGenerator::Pointer threadLocalInstance = [] {
    MersenneTwisterRandomVariateGenerator::Pointer instance = MersenneTwisterRandomVariateGenerator::CreateInstance();
    instance->SetSeed();
    return instance;
  }();

  return threadLocalInstance;

[N-Dekker]

The alternative is that each algorithm initializes its own random number generate ( either seeded by a user provided value, clock or the next seed from the global ). Using an algorithm specific generator removes the dependency on the global ( even a thread local is still a global and could cause issues is some cases, especially with certain thread pool or task based threading models.)

Interesting, @blowekamp . I guess for elastix we should then replace or rewrite ITK's random iterators, right? elastix's ImageRandomSampler currently still uses ITK's ImageRandomConstIteratorWithIndex, for example, which does MersenneTwisterRandomVariateGenerator::New() during its initialization. MersenneTwisterRandomVariateGenerator::New() calls GetNextSeed(), adjusting the global (shared) state. 🤷

[blowekamp]

The alternative is that each algorithm initializes its own random number generate ( either seeded by a user provided value, clock or the next seed from the global ). Using an algorithm specific generator removes the dependency on the global ( even a thread local is still a global and could cause issues is some cases, especially with certain thread pool or task based threading models.)

Interesting, @blowekamp . I guess for elastix we should then replace or rewrite ITK's random iterators, right? elastix's ImageRandomSampler currently still uses ITK's ImageRandomConstIteratorWithIndex, for example, which does MersenneTwisterRandomVariateGenerator::New() during its initialization. MersenneTwisterRandomVariateGenerator::New() calls GetNextSeed(), adjusting the global (shared) state. 🤷

So the initial issue is that if you run the algorithm (elastix) concurrently in multiple threads without setting the seed, then the behavior is randomized? Are there seeds available to set per algorithm or the only option is setting a global seed?

Have you looked ImageRegistrationMethodv4? A seed can be set for each instance of the class, and it uses it own MersenneTwisterRandomVariateGenerator to generate seeds for each of its components.

Additionally, you may be interested in the method use the the Image Noise filters e.g.(https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Filtering/ImageNoise/include/itkSaltAndPepperNoiseImageFilter.hxx#L46-L54). Here the filter has a seed, and each thread modifies the seed for different per-seed values. This has the effect of each thread having a different random seed, but is reproducible by not cause a non-deterministic race condition for the seed. However the results do very with the number of threads used.

[N-Dekker]

@blowekamp Even when each thread has its own instance of MersenneTwisterRandomVariateGenerator, they still share the one static variable m_StaticDiffer:

ITK/Modules/Core/Common/src/itkMersenneTwisterRandomVariateGenerator.cxx

Line 39 in db10423

std::atomic<MersenneTwisterRandomVariateGenerator::IntegerType> m_StaticDiffer{};

I guess that might still cause nondeterministic behavior, right? Even if it is "thread-safe" (atomic), it may still affect the behavior when one thread does m_StaticDiffer++ before or after another thread. 🤷

[blowekamp]

I guess that might still cause nondeterministic behavior, right? Even if it is "thread-safe" (atomic), it may still affect the behavior when one thread does m_StaticDiffer++ before or after another thread. 🤷

My understanding is that m_StaticDiffer is use to set the seeds sequentially. My recommendations above are to explicitly set the seeds for the algorithm, and not rely on the global state.

[N-Dekker]

My understanding is that m_StaticDiffer is use to set the seeds sequentially. My recommendations above are to explicitly set the seeds for the algorithm, and not rely on the global state.

Thanks @blowekamp, I see now, two different generators produce the same random number sequence, when they have the same initial seed:

using itk::Statistics::MersenneTwisterRandomVariateGenerator;

const auto generator1 = MersenneTwisterRandomVariateGenerator::New();
const auto generator2 = MersenneTwisterRandomVariateGenerator::New();

generator1->SetSeed(42);
generator2->SetSeed(42);

std::cout << generator1->GetIntegerVariate() << '\n'
          << generator2->GetIntegerVariate() << '\n'
          << generator1->GetIntegerVariate() << '\n'
          << generator2->GetIntegerVariate() << '\n'
          << generator1->GetIntegerVariate() << '\n'
          << generator2->GetIntegerVariate() << '\n';

Output:

1608637542
1608637542
3421126067
3421126067
4083286876
4083286876

👍

I think it's a bit unfortunate that the global m_StaticDiffer is still modified during the creation of those generators. Maybe it would be nice to have an extra MersenneTwisterRandomVariateGenerator::Create function that would not modify the global data. What do you think?

[blowekamp]

I think it's a bit unfortunate that the global m_StaticDiffer is still modified during the creation of those generators. Maybe it would be nice to have an extra MersenneTwisterRandomVariateGenerator::Create function that would not modify the global data. What do you think?

Yes. That makes sense to me.

[dzenanz]

From the discussion here, it looks like this PR is not a good solution. Have we identified an alternative solution? If so, who should propose a PR? Or should this PR be abandoned?

[N-Dekker]

Have we identified an alternative solution? If so, who should propose a PR? Or should this PR be abandoned?

I'm still looking into elastix, trying to replace all those GetInstance() calls, but it's not so easy 🤷 I will comment here again when it works (if it works).

[N-Dekker]

@michalmeszaroswork Thanks again for your work on this topic. Do you agree that with the recent merge of pull request SuperElastix/elastix#1323, your ITK pull request is no longer necessary, for elastix?

[hjmjohnson]

Fixes made in SuperElastix/elastix#1323 make this PR unnecessary.