EigenRandom

include/flucoma/algorithms/util/EigenRandom.hpp

@@ -0,0 +1,162 @@
+#include "../../data/FluidIndex.hpp"
+#include <Eigen/Core>
+#include <cmath>
+#include <optional>
+#include <random>
+
+namespace fluid::algorithm {
+
+struct RandomSeed
+{
+  RandomSeed(index seed) : mSeed{seed} {}
+  RandomSeed() : RandomSeed(-1) {}
+  operator std::optional<size_t>()
+  {
+    return mSeed >= 0 ? mSeed : std::optional<index>{};
+  }
+
+  index mSeed;
+};
+
+
+template <typename T>
+struct is_complex : std::false_type
+{};
+template <typename U>
+struct is_complex<std::complex<U>> : std::true_type
+{};
+template <typename T>
+inline constexpr bool is_complex_v = is_complex<T>::value;
+
+template <typename T>
+struct value_type_of
+{
+  using type = T;
+};
+
+template <typename T>
+struct value_type_of<std::complex<T>>
+{
+  using type = T;
+};
+
+template <typename T>
+using value_type_of_t = typename value_type_of<T>::type;
+
+template <typename T>
+struct Range
+{
+  T min;
+  T max;
+};
+
+// deduction guide
+template <typename T>
+Range(T, T) -> Range<T>;
+
+template <typename T>
+Range<T> DefaultRange =
+    Range<T>{std::numeric_limits<T>::min(), std::numeric_limits<T>::max()};
+
+template <>
+auto DefaultRange<double> = Range<double>{-1.0, 1.0};
+
+template <>
+auto DefaultRange<float> = Range<double>{-1.0f, 1.0f};
+
+
+// Wraps up C++ random stuff to suit our needs
+// can cope with uniform distributions of integral or real or complex real types
+// replciates Eigen::Random default ranges (-1:1 for reals, -max:max for
+// integral types) but also allows settable ranges (unlike Eigen::Random)
+template <typename T,
+          typename D = std::conditional_t<
+              std::is_integral_v<T>, std::uniform_int_distribution<T>,
+              std::uniform_real_distribution<value_type_of_t<T>>>>
+struct RandomGenerator
+{
+  RandomGenerator(std::optional<size_t> seed, Range<value_type_of_t<T>> range)
+      : g{seed ? *seed : rd()}, d{range.min, range.max}
+  {}
+  RandomGenerator(std::optional<size_t> seed)
+      : RandomGenerator(seed, DefaultRange<value_type_of_t<T>>)
+  {}
+  RandomGenerator(const RandomGenerator& x) : g{x.g}, d{x.d} {}
+
+  T operator()() const
+  {
+    // if complex, we need to generate two random values
+    if constexpr (!is_complex_v<T>)
+      return d(g);
+    else
+      return T{d(g), d(g)};
+  }
+
+private:
+  std::random_device      rd;
+  mutable std::mt19937_64 g;
+
+  mutable D d;
+};
+
+
+template <typename EigenType>
+auto EigenRandom(index rows, index cols, RandomSeed seed,
+                 Range<value_type_of_t<typename EigenType::Scalar>> range)
+{
+  return EigenType::NullaryExpr(
+      rows, cols, RandomGenerator<typename EigenType::Scalar>{seed, range});
+}
+
+template <typename EigenType>
+auto EigenRandom(index rows, index cols, RandomSeed seed)
+{
+  using T = value_type_of_t<typename EigenType::Scalar>;
+  return EigenRandom<EigenType>(rows, cols, seed, DefaultRange<T>);
+}
+
+template <typename EigenType>
+auto EigenRandom(index size, RandomSeed seed,
+                 Range<value_type_of_t<typename EigenType::Scalar>> range)
+{
+  static_assert(EigenType::RowsAtCompileTime == 1 ||
+                    EigenType::ColsAtCompileTime == 1,
+                "Eigen Type is not a vector, did you mean to call (rows,cols) "
+                "version?");
+  if constexpr (EigenType::ColsAtCompileTime == 1)
+  {
+    return EigenType::NullaryExpr(
+        size, 1, RandomGenerator<typename EigenType::Scalar>{seed, range});
+  }
+  if constexpr (EigenType::RowsAtCompileTime == 1)
+  {
+    return EigenType::NullaryExpr(
+        1, size, RandomGenerator<typename EigenType::Scalar>{seed, range});
+  }
+}
+
+template <typename EigenType>
+auto EigenRandom(index size, RandomSeed seed)
+{
+  using T = value_type_of_t<typename EigenType::Scalar>;
+  return EigenRandom<EigenType>(size, seed, DefaultRange<T>);
+}
+
+template <typename EigenType>
+auto EigenRandomPhase(index rows, index cols, RandomSeed seed)
+{
+
+  static_assert(is_complex_v<typename EigenType::Scalar>,
+                "This only works for complex valued Eigen types");
+
+  constexpr double twopi = 2 * M_PI;
+  using T = value_type_of_t<typename EigenType::Scalar>;
+  auto f = [rg = RandomGenerator<T>{seed, Range{0.0, twopi}}]() {
+    return std::polar(1.0, rg());
+  };
+
+  return EigenType::NullaryExpr(rows, cols, f);
+}
+
+
+} // namespace fluid::algorithm

tests/CMakeLists.txt

@@ -123,6 +123,9 @@ target_link_libraries(TestEnvelopeSeg PRIVATE TestSignals)
 target_link_libraries(TestEnvelopeGate PRIVATE TestSignals)
 target_link_libraries(TestTransientSlice PRIVATE TestSignals)
 
+add_test_executable(TestEigenRandom algorithms/util/TestEigenRandom.cpp)
+
+
 include(CTest)
 include(Catch)
 
@@ -149,5 +152,6 @@ catch_discover_tests(TestBufferedProcess WORKING_DIRECTORY "${CMAKE_BINARY_DIR}"
 
 catch_discover_tests(TestMLP WORKING_DIRECTORY "${CMAKE_BINARY_DIR}")
 catch_discover_tests(TestKMeans WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
+catch_discover_tests(TestEigenRandom WORKING_DIRECTORY ${CMAKE_BINARY_DIR})
 
 add_compile_tests("FluidTensor Compilation Tests" data/compile_tests/TestFluidTensor_Compile.cpp) 

tests/algorithms/util/TestEigenRandom.cpp

@@ -0,0 +1,127 @@
+#define CATCH_CONFIG_MAIN
+#include <Eigen/Core>
+#include <catch2/catch_all.hpp>
+#include <flucoma/algorithms/util/EigenRandom.hpp>
+#include <complex>
+#include <iostream>
+#include <optional>
+#include <random>
+
+namespace fluid::algorithm {
+
+TEST_CASE("EigenRandom Can Make Basic Random Containers")
+{
+
+  SECTION("Matrix of Double")
+  {
+    Eigen::MatrixXd mrnd = EigenRandom<Eigen::MatrixXd>(3, 4, RandomSeed{});
+
+    // Right size
+    REQUIRE((mrnd.rows() == 3 && mrnd.cols() == 4));
+    // No NaNs or infs
+    REQUIRE((mrnd.array().isFinite()).all());
+    // Not constant value
+    auto first = mrnd(0, 0);
+    REQUIRE_FALSE((mrnd.array() == first).all());
+  }
+
+  SECTION("2D Array of Double")
+  {
+    Eigen::ArrayXXd mrnd = EigenRandom<Eigen::ArrayXXd>(3, 4, RandomSeed{});
+
+    // Right size
+    REQUIRE((mrnd.rows() == 3 && mrnd.cols() == 4));
+    // No NaNs or infs
+    REQUIRE((mrnd.array().isFinite()).all());
+    // Not constant value
+    auto first = mrnd(0, 0);
+    REQUIRE_FALSE((mrnd.array() == first).all());
+  }
+
+  SECTION("Vector of Double")
+  {
+    Eigen::VectorXd vrnd = EigenRandom<Eigen::VectorXd>(4, RandomSeed{});
+
+    // Right size
+    REQUIRE((vrnd.size() == 4));
+    // No NaNs or infs
+    REQUIRE((vrnd.array().isFinite()).all());
+    // Not constant value
+    auto first = vrnd(0);
+    REQUIRE_FALSE((vrnd.array() == first).all());
+  }
+
+  SECTION("Vector of int")
+  {
+    Eigen::VectorXi vrnd = EigenRandom<Eigen::VectorXi>(4, RandomSeed{});
+
+    // Right size
+    REQUIRE((vrnd.size() == 4));
+    // Not constant value
+    auto first = vrnd(0);
+    REQUIRE_FALSE((vrnd.array() == first).all());
+  }
+
+  SECTION("1D Array of Double")
+  {
+    Eigen::ArrayXd mrnd = EigenRandom<Eigen::ArrayXd>(4, RandomSeed{});
+
+    // Right size
+    REQUIRE(mrnd.size() == 4);
+    // No NaNs or infs
+    REQUIRE((mrnd.array().isFinite()).all());
+    // Not constant value
+    auto first = mrnd(0);
+    REQUIRE_FALSE((mrnd.array() == first).all());
+  }
+
+  SECTION("2D Array of Complex Double")
+  {
+    Eigen::ArrayXXcd mrnd = EigenRandom<Eigen::ArrayXXcd>(3, 4, RandomSeed{});
+
+    // Right size
+    REQUIRE((mrnd.rows() == 3 && mrnd.cols() == 4));
+    // No NaNs or infs
+    REQUIRE((mrnd.array().isFinite()).all());
+    // Not constant value
+    auto first = mrnd(0, 0);
+    REQUIRE_FALSE((mrnd.array() == first).all());
+  }
+}
+
+TEST_CASE("EigenRandom Can Manually Seet Random Seed For Repeatability")
+{
+
+  Eigen::MatrixXd mrnd1 = EigenRandom<Eigen::MatrixXd>(3, 4, RandomSeed{42});
+  Eigen::MatrixXd mrnd2 = EigenRandom<Eigen::MatrixXd>(3, 4, RandomSeed{42});
+
+  // Same seed -> same result
+  REQUIRE((mrnd1.array() == mrnd2.array()).all());
+
+  // Different seed -> different result
+  Eigen::MatrixXd mrnd3 = EigenRandom<Eigen::MatrixXd>(3, 4, RandomSeed{4203});
+  REQUIRE_FALSE((mrnd1.array() == mrnd3.array()).all());
+}
+
+TEST_CASE("Eigen Random Can Set Min and Max")
+{
+  double          min = -0.4;
+  double          max = 0.1;
+  Eigen::MatrixXd mrnd =
+      EigenRandom<Eigen::MatrixXd>(3, 4, RandomSeed{42}, Range{-0.4, 0.1});
+
+  REQUIRE((mrnd.minCoeff() >= min && mrnd.maxCoeff() < max));
+}
+
+TEST_CASE("EigenRandom Can make complex array of random phase, unit magnitude")
+{
+  Eigen::ArrayXXcd a = EigenRandomPhase<Eigen::ArrayXXcd>(4, 3, RandomSeed{});
+  constexpr double twopi = 2 * M_PI;
+
+  std::cout << a.arg();
+  REQUIRE((a.arg().maxCoeff() - a.arg().minCoeff() <= twopi));
+  REQUIRE_THAT(a.abs().reshaped(),
+               Catch::Matchers::AllMatch(Catch::Matchers::WithinRel(1.0)));
+}
+
+} // namespace fluid::algorithm