ArrayHandleRandomStandardNormal.h

Go to the documentation of this file.

//============================================================================
//  The contents of this file are covered by the Viskores license. See
//  LICENSE.txt for details.
//
//  By contributing to this file, all contributors agree to the Developer
//  Certificate of Origin Version 1.1 (DCO 1.1) as stated in DCO.txt.
//============================================================================
 
//============================================================================
//  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 viskores_count_ArrayHandleRandomStandardNormal_h
#define viskores_count_ArrayHandleRandomStandardNormal_h
 
#include <viskores/Math.h>
#include <viskores/cont/ArrayHandleRandomUniformReal.h>
#include <viskores/cont/ArrayHandleTransform.h>
#include <viskores/cont/ArrayHandleZip.h>
 
namespace viskores
{
namespace cont
{
namespace detail
{
struct BoxMuller
{
  VISKORES_EXEC_CONT viskores::Float32 operator()(
    const viskores::Pair<viskores::Float32, viskores::Float32>& uv) const
  {
    // We take two U(0, 1) and return one N(0, 1)
    return viskores::Sqrt(-2.0f * viskores::Log(uv.first)) *
      viskores::Cos(2.0f * viskores::TwoPif() * uv.second);
  }
 
  VISKORES_EXEC_CONT viskores::Float64 operator()(
    const viskores::Pair<viskores::Float64, viskores::Float64>& uv) const
  {
    // We take two U(0, 1) and return one N(0, 1)
    return viskores::Sqrt(-2.0 * viskores::Log(uv.first)) *
      viskores::Cos(2 * viskores::TwoPi() * uv.second);
  }
};
} //detail
 
template <typename Real = viskores::Float64>
class VISKORES_ALWAYS_EXPORT ArrayHandleRandomStandardNormal
  : public viskores::cont::ArrayHandleTransform<
      viskores::cont::ArrayHandleZip<viskores::cont::ArrayHandleRandomUniformReal<Real>,
                                     viskores::cont::ArrayHandleRandomUniformReal<Real>>,
      detail::BoxMuller>
{
public:
  using SeedType = viskores::Vec<viskores::UInt32, 1>;
  using UniformReal = viskores::cont::ArrayHandleRandomUniformReal<Real>;
 
  VISKORES_ARRAY_HANDLE_SUBCLASS(
    ArrayHandleRandomStandardNormal,
    (ArrayHandleRandomStandardNormal<Real>),
    (viskores::cont::ArrayHandleTransform<
      viskores::cont::ArrayHandleZip<viskores::cont::ArrayHandleRandomUniformReal<Real>,
                                     viskores::cont::ArrayHandleRandomUniformReal<Real>>,
      detail::BoxMuller>));
 
  explicit ArrayHandleRandomStandardNormal(viskores::Id length,
                                           SeedType seed = { std::random_device{}() })
    : Superclass(viskores::cont::make_ArrayHandleZip(UniformReal{ length, seed },
                                                     UniformReal{ length, { ~seed[0] } }),
                 detail::BoxMuller{})
  {
  }
};
}
}
#endif // viskores_count_ArrayHandleRandomStandardNormal_h