Merge topic 'philox'

ef49a71aa Merge branch 'philox' of gitlab.kitware.com:ollielo/vtk-m into philox
563642400 remove outdated comments on type-punning
be7c18544 tried to supress warning about type punning
c795f74b2 Merge branch 'master' into philox
2b43fcb25 change function signature for work with MSVC
13e7ac362 Merge branch 'philox' of gitlab.kitware.com:ollielo/vtk-m into philox
1e2fd540d add Doxygen documentation
f2ea17917 add Doxygen documentation
...

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !1990

vtkm/CMakeLists.txt

@@ -113,3 +113,6 @@ add_subdirectory(io)
 
 #add the source folder
 add_subdirectory(source)
+
+#add Pseudo Random Number Generator folder
+add_subdirectory(random)

vtkm/cont/ArrayHandleRandomUniformBits.h

@@ -0,0 +1,98 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+#ifndef vtk_m_cont_ArrayHandlePhiloxURBG_h
+#define vtk_m_cont_ArrayHandlePhiloxURBG_h
+
+#include <random>
+#include <vtkm/cont/ArrayHandleImplicit.h>
+#include <vtkm/random/Philox.h>
+
+namespace vtkm
+{
+namespace cont
+{
+
+namespace detail
+{
+struct PhiloxFunctor
+{
+  using SeedType = vtkm::Vec<vtkm::UInt32, 1>;
+
+  PhiloxFunctor() = default;
+
+  explicit PhiloxFunctor(SeedType seed)
+    : Seed(seed)
+  {
+  }
+
+  VTKM_EXEC_CONT
+  vtkm::UInt64 operator()(vtkm::Id index) const
+  {
+    using philox_functor = vtkm::random::PhiloxFunctor2x32x10;
+    using counters_type = typename philox_functor::counters_type;
+
+    auto idx = static_cast<vtkm::UInt64>(index);
+    counters_type counters{ static_cast<vtkm::UInt32>(idx), static_cast<vtkm::UInt32>(idx >> 32) };
+    counters_type result = philox_functor{}(counters, Seed);
+    return static_cast<vtkm::UInt64>(result[0]) | static_cast<vtkm::UInt64>(result[1]) << 32;
+  }
+
+private:
+  const SeedType Seed{};
+}; // class PhiloxFunctor
+} // namespace detail
+
+/// \brief An \c ArrayHandle that provides a source of random bits
+///
+/// \c ArrayHandleRandomUniformBits is a specialization of ArrayHandleImplicit.
+/// It takes a user supplied seed and hash it with the a given index value. The
+/// hashed value is the value of the array at that position.
+///
+/// Currently, Philox2x32x10 as described in the
+///   "Parallel Random Numbers: As Easy as 1, 2, 3," Proceedings of the
+///   International Conference for High Performance Computing, Networking,
+///   Storage and Analysis (SC11)
+/// is used as the hash function.
+///
+/// Note: In contrast to traditional random number generator,
+/// ArrayHandleRandomUniformBits does not have "state", i.e. multiple calls
+/// the Get() method with the same index will always return the same hash value.
+/// To ge a new set of random bits, create a new ArrayHandleRandomUniformBits
+/// with a different seed.
+class VTKM_ALWAYS_EXPORT ArrayHandleRandomUniformBits
+  : public vtkm::cont::ArrayHandleImplicit<detail::PhiloxFunctor>
+{
+public:
+  using SeedType = vtkm::Vec<vtkm::UInt32, 1>;
+
+  VTKM_ARRAY_HANDLE_SUBCLASS_NT(ArrayHandleRandomUniformBits,
+                                (vtkm::cont::ArrayHandleImplicit<detail::PhiloxFunctor>));
+
+  /// The type of seed is specifically designed to be an vtkm::Vec<> to provide
+  /// type safety for the parameters so user will not transpose two integer parameters.
+  explicit ArrayHandleRandomUniformBits(vtkm::Id length, SeedType seed = { std::random_device{}() })
+    : Superclass(detail::PhiloxFunctor(seed), length)
+  {
+  }
+}; // class ArrayHandleRandomUniformBits
+}
+} // namespace vtkm::cont
+
+//=============================================================================
+// Specializations of serialization related classes
+/// @cond SERIALIZATION
+
+namespace vtkm
+{
+namespace cont
+{
+}
+} // namespace vtkm::cont
+#endif //vtk_m_cont_ArrayHandlePhiloxURBG_h

vtkm/cont/CMakeLists.txt

@@ -32,6 +32,7 @@ set(headers
   ArrayHandleMultiplexer.h
   ArrayHandlePermutation.h
   ArrayHandleReverse.h
+  ArrayHandleRandomUniformBits.h
   ArrayHandleSOA.h
   ArrayHandleSwizzle.h
   ArrayHandleTransform.h

vtkm/cont/testing/CMakeLists.txt

@@ -47,6 +47,7 @@ set(unit_tests
   UnitTestArrayHandleIndex.cxx
   UnitTestArrayHandleReverse.cxx
   UnitTestArrayHandlePermutation.cxx
+  UnitTestArrayHandleRandomUniformBits.cxx
   UnitTestArrayHandleSwizzle.cxx
   UnitTestArrayHandleThreadSafety.cxx
   UnitTestArrayHandleTransform.cxx

vtkm/cont/testing/UnitTestArrayHandleRandomUniformBits.cxx

@@ -0,0 +1,44 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+
+#include <vtkm/cont/ArrayHandleRandomUniformBits.h>
+#include <vtkm/cont/testing/Testing.h>
+
+void TestArrayHandleRandomUniformBits()
+{
+  auto actual0 = vtkm::cont::ArrayHandleRandomUniformBits(10, { 0 });
+  // result from Random123 sample implementation of philox2x32x10
+  std::vector<vtkm::UInt64> expected0{ 0x6cd10df2ff1dae59, 0x5f3adb6bdcdce855, 0x3fbb6394049f6998,
+                                       0xbd592d1202a74512, 0x8a115b62c08084ef, 0x1411803b3bb7eefa,
+                                       0x7d138a2280027d0e, 0x318a7703a1da82c5, 0xdcd79c6998975579,
+                                       0x6cb1a07c91f81109 };
+
+
+  auto result =
+    vtkm::cont::testing::test_equal_ArrayHandles(actual0, vtkm::cont::make_ArrayHandle(expected0));
+  VTKM_TEST_ASSERT(result, result.GetMergedMessage());
+
+  // initialize with seed = 100, could be "iteration number" in actual use case.
+  auto actual100 = vtkm::cont::ArrayHandleRandomUniformBits(10, { 100 });
+  // result from Random123 sample implementation of philox2x32x10
+  std::vector<vtkm::UInt64> expected100{
+    0xbd35360836122ea3, 0xe033b74acce7aa5f, 0xc0fbb65cba93ecd7, 0xe3fee2812b77e480,
+    0x92e5c7d563767971, 0xd99e952fb054fc19, 0xb8f2adc12094ad29, 0xb7dcb35fea8c27ac,
+    0x9c7b779e88270c45, 0x7325b123dc32e01d,
+  };
+  auto result100 = vtkm::cont::testing::test_equal_ArrayHandles(
+    actual100, vtkm::cont::make_ArrayHandle(expected100));
+  VTKM_TEST_ASSERT(result, result.GetMergedMessage());
+}
+
+int UnitTestArrayHandleRandomUniformBits(int argc, char* argv[])
+{
+  return vtkm::cont::testing::Testing::Run(TestArrayHandleRandomUniformBits, argc, argv);
+}

vtkm/random/CMakeLists.txt

@@ -0,0 +1,17 @@
+##============================================================================
+##  Copyright (c) Kitware, Inc.
+##  All rights reserved.
+##  See LICENSE.txt for details.
+##
+##  This software is distributed WITHOUT ANY WARRANTY; without even
+##  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+##  PURPOSE.  See the above copyright notice for more information.
+##============================================================================
+
+set(headers
+  Philox.h)
+
+vtkm_declare_headers(${headers})
+
+#-----------------------------------------------------------------------------
+add_subdirectory(testing)

vtkm/random/Philox.h

@@ -0,0 +1,124 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+#ifndef vtk_m_random_Philox_h
+#define vtk_m_random_Philox_h
+
+#include <vtkm/Types.h>
+
+namespace vtkm
+{
+namespace random
+{
+namespace detail
+{
+VTKM_EXEC_CONT vtkm::Vec<vtkm::UInt32, 2> mulhilo(vtkm::UInt32 a, vtkm::UInt32 b)
+{
+  vtkm::UInt64 r = static_cast<vtkm::UInt64>(a) * b;
+  auto lo = static_cast<vtkm::UInt32>(r);
+  vtkm::UInt32 hi = static_cast<vtkm::UInt32>(r >> 32);
+  return { lo, hi };
+}
+
+template <typename UIntType, std::size_t N, UIntType... consts>
+struct philox_parameters;
+
+template <typename T, T M0, T C0>
+struct philox_parameters<T, 2, M0, C0>
+{
+  static constexpr Vec<T, 1> multipliers = { M0 };
+  static constexpr Vec<T, 1> round_consts = { C0 };
+};
+
+template <typename T, T M0, T C0, T M1, T C1>
+struct philox_parameters<T, 4, M0, C0, M1, C1>
+{
+  static constexpr vtkm::Vec<T, 2> multipliers = { M0, M1 };
+  static constexpr vtkm::Vec<T, 2> round_consts = { C0, C1 };
+};
+
+template <typename UIntType, std::size_t N, std::size_t R, UIntType... consts>
+class philox_functor;
+
+template <typename UIntType, std::size_t R, UIntType... consts>
+class philox_functor<UIntType, 2, R, consts...>
+{
+public:
+  using counters_type = vtkm::Vec<UIntType, 2>;
+  using keys_type = vtkm::Vec<UIntType, 1>;
+
+  VTKM_EXEC_CONT counters_type operator()(counters_type counters, keys_type keys) const
+  {
+    for (std::size_t i = 0; i < R; ++i)
+    {
+      counters = round(counters, keys);
+      keys = bump_keys(keys);
+    }
+    return counters;
+  }
+
+private:
+  static VTKM_EXEC_CONT counters_type round(counters_type counters, keys_type round_keys)
+  {
+    auto constexpr multipliers = philox_parameters<UIntType, 2, consts...>::multipliers;
+    vtkm::Vec<UIntType, 2> r = mulhilo(multipliers[0], counters[0]);
+    return { r[1] ^ round_keys[0] ^ counters[1], r[0] };
+  }
+
+  static VTKM_EXEC_CONT keys_type bump_keys(keys_type keys)
+  {
+    auto constexpr round_consts = philox_parameters<UIntType, 2, consts...>::round_consts;
+    return { keys[0] + round_consts[0] };
+  }
+};
+
+template <typename UIntType, std::size_t R, UIntType... consts>
+class philox_functor<UIntType, 4, R, consts...>
+{
+  using counters_type = vtkm::Vec<UIntType, 4>;
+  using keys_type = vtkm::Vec<UIntType, 2>;
+
+  static VTKM_EXEC_CONT counters_type round(counters_type counters, keys_type round_keys)
+  {
+    auto constexpr multipliers = philox_parameters<UIntType, 4, consts...>::multipliers;
+    vtkm::Vec<UIntType, 2> r0 = mulhilo(multipliers[0], counters[0]);
+    vtkm::Vec<UIntType, 2> r1 = mulhilo(multipliers[1], counters[2]);
+    return {
+      r1[1] ^ round_keys[0] ^ counters[1], r1[0], r0[1] ^ round_keys[1] ^ counters[3], r0[0]
+    };
+  }
+
+  static VTKM_EXEC_CONT keys_type bump_key(keys_type keys)
+  {
+    auto constexpr round_consts = philox_parameters<UIntType, 4, consts...>::round_consts;
+    keys[0] += round_consts[0];
+    keys[1] += round_consts[1];
+    return keys;
+  }
+
+public:
+  VTKM_EXEC_CONT counters_type operator()(counters_type counters, keys_type keys) const
+  {
+    for (std::size_t i = 0; i < R; ++i)
+    {
+      counters = round(counters, keys);
+      keys = bump_key(keys);
+    }
+    return counters;
+  }
+};
+
+} // namespace detail
+
+using PhiloxFunctor2x32x7 = detail::philox_functor<vtkm::UInt32, 2, 7, 0xD256D193, 0x9E3779B9>;
+using PhiloxFunctor2x32x10 = detail::philox_functor<vtkm::UInt32, 2, 10, 0xD256D193, 0x9E3779B9>;
+
+} // namespace random
+} // namespace vtkm
+#endif //vtk_m_random_Philox_h

vtkm/random/testing/CMakeLists.txt

@@ -0,0 +1,15 @@
+##============================================================================
+##  Copyright (c) Kitware, Inc.
+##  All rights reserved.
+##  See LICENSE.txt for details.
+##
+##  This software is distributed WITHOUT ANY WARRANTY; without even
+##  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+##  PURPOSE.  See the above copyright notice for more information.
+##============================================================================
+
+set(unit_tests
+  UnitTestPhiloxRNG.cxx)
+
+vtkm_unit_tests(SOURCES ${unit_tests}
+  ALL_BACKENDS)

vtkm/random/testing/UnitTestPhiloxRNG.cxx

@@ -0,0 +1,49 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+
+#include <vtkm/cont/testing/Testing.h>
+#include <vtkm/random/Philox.h>
+
+void TestPhiloxRNG2x32x7()
+{
+  vtkm::random::PhiloxFunctor2x32x7 f2x32x7;
+  using counters_type = typename vtkm::random::PhiloxFunctor2x32x7::counters_type;
+
+  // test cases from original Random123 implementation
+  VTKM_TEST_ASSERT(f2x32x7({}, {}) == counters_type{ 0x257a3673, 0xcd26be2a });
+  VTKM_TEST_ASSERT(f2x32x7({ 0xffffffff, 0xffffffff }, { 0xffffffff }) ==
+                   counters_type{ 0xab302c4d, 0x3dc9d239 });
+  VTKM_TEST_ASSERT(f2x32x7({ 0x243f6a88, 0x85a308d3 }, { 0x13198a2e }) ==
+                   counters_type{ 0xbedbbe6b, 0xe4c770b3 });
+}
+
+void TestPhiloxRNG2x32x10()
+{
+  vtkm::random::PhiloxFunctor2x32x10 f2x32x10;
+  using counters_type = typename vtkm::random::PhiloxFunctor2x32x10::counters_type;
+
+  // test cases from original Random123 implementation
+  VTKM_TEST_ASSERT(f2x32x10({}, {}) == counters_type{ 0xff1dae59, 0x6cd10df2 });
+  VTKM_TEST_ASSERT(f2x32x10({ 0xffffffff, 0xffffffff }, { 0xffffffff }) ==
+                   counters_type{ 0x2c3f628b, 0xab4fd7ad });
+  VTKM_TEST_ASSERT(f2x32x10({ 0x243f6a88, 0x85a308d3 }, { 0x13198a2e }) ==
+                   counters_type{ 0xdd7ce038, 0xf62a4c12 });
+}
+
+void TestPhiloxRNG()
+{
+  TestPhiloxRNG2x32x7();
+  TestPhiloxRNG2x32x10();
+}
+
+int UnitTestPhiloxRNG(int argc, char* argv[])
+{
+  return vtkm::cont::testing::Testing::Run(TestPhiloxRNG, argc, argv);
+}