diff --git a/vtkm/CMakeLists.txt b/vtkm/CMakeLists.txt index c73635cc8..e6d54c639 100644 --- a/vtkm/CMakeLists.txt +++ b/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) diff --git a/vtkm/cont/ArrayHandleRandomUniformBits.h b/vtkm/cont/ArrayHandleRandomUniformBits.h new file mode 100644 index 000000000..eadf4500f --- /dev/null +++ b/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 +#include +#include + +namespace vtkm +{ +namespace cont +{ + +namespace detail +{ +struct PhiloxFunctor +{ + using SeedType = vtkm::Vec; + + 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(index); + counters_type counters{ static_cast(idx), static_cast(idx >> 32) }; + counters_type result = philox_functor{}(counters, Seed); + return static_cast(result[0]) | static_cast(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 +{ +public: + using SeedType = vtkm::Vec; + + VTKM_ARRAY_HANDLE_SUBCLASS_NT(ArrayHandleRandomUniformBits, + (vtkm::cont::ArrayHandleImplicit)); + + /// 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 diff --git a/vtkm/cont/CMakeLists.txt b/vtkm/cont/CMakeLists.txt index 980ab4172..3cf9dd1b6 100644 --- a/vtkm/cont/CMakeLists.txt +++ b/vtkm/cont/CMakeLists.txt @@ -32,6 +32,7 @@ set(headers ArrayHandleMultiplexer.h ArrayHandlePermutation.h ArrayHandleReverse.h + ArrayHandleRandomUniformBits.h ArrayHandleSOA.h ArrayHandleSwizzle.h ArrayHandleTransform.h diff --git a/vtkm/cont/testing/CMakeLists.txt b/vtkm/cont/testing/CMakeLists.txt index 45ae558d4..0fb27dc53 100644 --- a/vtkm/cont/testing/CMakeLists.txt +++ b/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 diff --git a/vtkm/cont/testing/UnitTestArrayHandleRandomUniformBits.cxx b/vtkm/cont/testing/UnitTestArrayHandleRandomUniformBits.cxx new file mode 100644 index 000000000..2374b0fc1 --- /dev/null +++ b/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 +#include + +void TestArrayHandleRandomUniformBits() +{ + auto actual0 = vtkm::cont::ArrayHandleRandomUniformBits(10, { 0 }); + // result from Random123 sample implementation of philox2x32x10 + std::vector 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 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); +} diff --git a/vtkm/random/CMakeLists.txt b/vtkm/random/CMakeLists.txt new file mode 100644 index 000000000..9ab10ff01 --- /dev/null +++ b/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) diff --git a/vtkm/random/Philox.h b/vtkm/random/Philox.h new file mode 100644 index 000000000..413d1933b --- /dev/null +++ b/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 + +namespace vtkm +{ +namespace random +{ +namespace detail +{ +VTKM_EXEC_CONT vtkm::Vec mulhilo(vtkm::UInt32 a, vtkm::UInt32 b) +{ + vtkm::UInt64 r = static_cast(a) * b; + auto lo = static_cast(r); + vtkm::UInt32 hi = static_cast(r >> 32); + return { lo, hi }; +} + +template +struct philox_parameters; + +template +struct philox_parameters +{ + static constexpr Vec multipliers = { M0 }; + static constexpr Vec round_consts = { C0 }; +}; + +template +struct philox_parameters +{ + static constexpr vtkm::Vec multipliers = { M0, M1 }; + static constexpr vtkm::Vec round_consts = { C0, C1 }; +}; + +template +class philox_functor; + +template +class philox_functor +{ +public: + using counters_type = vtkm::Vec; + using keys_type = vtkm::Vec; + + 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::multipliers; + vtkm::Vec 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::round_consts; + return { keys[0] + round_consts[0] }; + } +}; + +template +class philox_functor +{ + using counters_type = vtkm::Vec; + using keys_type = vtkm::Vec; + + static VTKM_EXEC_CONT counters_type round(counters_type counters, keys_type round_keys) + { + auto constexpr multipliers = philox_parameters::multipliers; + vtkm::Vec r0 = mulhilo(multipliers[0], counters[0]); + vtkm::Vec 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::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; +using PhiloxFunctor2x32x10 = detail::philox_functor; + +} // namespace random +} // namespace vtkm +#endif //vtk_m_random_Philox_h diff --git a/vtkm/random/testing/CMakeLists.txt b/vtkm/random/testing/CMakeLists.txt new file mode 100644 index 000000000..9158d0372 --- /dev/null +++ b/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) diff --git a/vtkm/random/testing/UnitTestPhiloxRNG.cxx b/vtkm/random/testing/UnitTestPhiloxRNG.cxx new file mode 100644 index 000000000..de7e40748 --- /dev/null +++ b/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 +#include + +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); +}