ochsnerd · ochsnerd · Jan 26, 2025 · Jan 4, 2025
diff --git a/Makefile b/Makefile
@@ -10,7 +10,7 @@ build:
 
 install: build interface.py
 	mkdir -p $(INSTALL_DIR)/bin
-	cp -n result/lerw $(INSTALL_DIR)/bin
+	cp result/lerw $(INSTALL_DIR)/bin
 	cp -n interface.py $(INSTALL_DIR)
 
 build_manual:

diff --git a/include/directions.hpp b/include/directions.hpp
@@ -3,6 +3,7 @@
 #include <algorithm>
 #include <array>
 #include <cstddef>
+#include <cstdint>
 #include <numeric>
 #include <random>
 
@@ -31,6 +32,95 @@ template <point Point> struct L2Direction {
   }
 };
 
+template <std::uniform_random_bit_generator RNG>
+constexpr auto random_sign(RNG &rng) -> std::int8_t {
+  using i = std::int8_t;
+  auto d = std::uniform_int_distribution<i>{0, 1};
+  return static_cast<i>(2 * d(rng) - 1);
+}
+
+// uniformly pick direction from surface of L1-ball
+template <point Point> struct L1Direction {
+  // find an integer solution to
+  // sum |x_i| = r
+  // by using the 'stars and bars' method.
+  // Example:
+  // For r=5, d=3
+  // |****|*
+  // corresponds to (0, 4, 1)
+  // After this, the sign of each coordinate is randomized.
+  // Because (-1 * 0 = 1 * 0 = 0), this oversamples points
+  // by a factor of 2 for every 0-coordinate. To combat
+  // this, reject a point with a probability of 0.5 for
+  // every 0-coordinate.
+  // TODO: The stars-and-bars could be adapted to get around this multi-counting
+
+  using result_type = Point;
+
+  using int_t = field<Point>::type;
+
+  static constexpr std::size_t d = dim<Point>();
+  static_assert(d > 0, "Out-of-bounds stuff happens when d = 0.");
+
+  template <std::uniform_random_bit_generator RNG>
+  constexpr auto operator()(int_t r, RNG &rng) -> Point {
+    auto proposed = get_integer_solution(r, rng);
+    while (reject_for_zeros(proposed.cbegin(), proposed.cend(), rng)) {
+      proposed = get_integer_solution(r, rng);
+    }
+    std::transform(proposed.cbegin(), proposed.cend(), proposed.begin(),
+                   [&rng](auto i) { return random_sign(rng) * i; });
+    return constructor<Point>{}(proposed.cbegin(), proposed.cend());
+  }
+
+  template <std::uniform_random_bit_generator RNG>
+  constexpr static auto get_integer_solution(int_t r,
+                                      RNG &rng) -> std::array<int_t, d> {
+    auto bars = sample_iota(r, rng);
+    bars[d - 1] = r + d;
+    std::sort(bars.begin(), bars.end());
+    auto stars = std::array<int_t, d>{};
+    std::adjacent_difference(bars.cbegin(), bars.cend(), stars.begin());
+    std::transform(stars.cbegin(), stars.cend(), stars.begin(),
+                   [&rng](auto s) { return (s - 1); });
+    return stars;
+  }
+
+  template <class InputIt, std::uniform_random_bit_generator RNG>
+  constexpr static auto reject_for_zeros(InputIt first, InputIt last,
+                                  RNG &rng) -> bool {
+    auto d = std::uniform_int_distribution<std::int8_t>{0, 1};
+    return std::transform_reduce(
+        first, last, false, std::logical_or{},
+        [&d, &rng](auto c) { return c == 0 && d(rng) == 0; });
+  }
+
+  template <std::uniform_random_bit_generator RNG>
+  constexpr static auto sample_iota(int_t r, RNG &rng) -> std::array<std::size_t, d> {
+    // equivalent to
+    // auto indices = std::vector<std::size_t>(r + d - 1);
+    // std::iota(indices.begin(), indices.end(), 1);
+    // auto result = std::array<std::size_t, d>{};
+    // std::sample(indices.cbegin(), indices.cend(), result.begin(), d - 1,
+    // rng);
+    // return result;
+    // but does this without materializing the indices-vector.
+    // This works fine for smallish d.
+    // See  Algorithm 3.4.2S of Seminumeric Algorithms (Knuth) for a less
+    // quadratic solution
+    auto dist = std::uniform_int_distribution<std::size_t>{1, r + d - 1};
+    auto result = std::array<std::size_t, d>{};
+    for (auto it = result.begin(); it + 1 < result.end(); ++it) {
+      auto next = dist(rng);
+      while (std::find(result.begin(), it, next) != it) {
+        next = dist(rng);
+      }
+      *it = next;
+    }
+    return result;
+  }
+};
+
 // uniformly pick direction from surface of linf-ball (cube centered at the
 // origin)
 template <point Point> struct LinfDirection {
@@ -53,26 +143,28 @@ template <point Point> struct LinfDirection {
 
   template <std::uniform_random_bit_generator RNG>
   constexpr auto operator()(int_t r, RNG &rng) -> Point {
-    r = std::abs(r);
     const auto k = choose_k(r, rng);
 
     auto coordinates = std::array<int_t, d>{};
     auto last = std::transform(
         coordinates.cbegin(), coordinates.cbegin() + k, coordinates.begin(),
-        [r, &rng](auto) mutable { return r * random_sign(rng); });
+        [r, &rng](auto) { return r * random_sign(rng); });
     constexpr int_t one{1}; // stupid typecasting
     auto dist = std::uniform_int_distribution<int_t>{one - r, r - one};
     std::transform(last, coordinates.end(), last,
-                   [&dist, &rng](auto) mutable { return dist(rng); });
+                   [&dist, &rng](auto) { return dist(rng); });
     std::shuffle(coordinates.begin(), coordinates.end(), rng);
     return constructor<Point>{}(coordinates.cbegin(), coordinates.cend());
   }
 
-  constexpr static auto A_k(std::uint16_t k, int_t r) -> float_t {
-    // boost forces floating point value return type, keep it throughout since
-    // for large r the number of points on faces dominate all others
-    const auto combinations = boost::math::binomial_coefficient<float_t>(d, k);
-    return combinations * std::pow(2, k) * std::pow(2 * r - 1, d - k);
+  template <std::uniform_random_bit_generator RNG>
+  constexpr static auto choose_k(int_t r, RNG &rng) -> std::uint16_t {
+    const auto ak_sums = A_k_sums(r);
+    const auto i =
+        std::uniform_real_distribution<float_t>{0, ak_sums[d - 1]}(rng);
+    return static_cast<std::uint16_t>(
+        std::distance(ak_sums.cbegin(),
+                      std::lower_bound(ak_sums.cbegin(), ak_sums.cend(), i)));
   }
 
   constexpr static auto A_k_sums(int_t r) -> std::array<float_t, d> {
@@ -88,20 +180,14 @@ template <point Point> struct LinfDirection {
     return ak_sums;
   }
 
-  template <std::uniform_random_bit_generator RNG>
-  constexpr static auto choose_k(int_t r, RNG &rng) -> std::uint16_t {
-    const auto ak_sums = A_k_sums(r);
-    const auto i =
-        std::uniform_real_distribution<float_t>{0, ak_sums[d - 1]}(rng);
-    return static_cast<std::uint16_t>(
-        std::distance(ak_sums.cbegin(),
-                      std::lower_bound(ak_sums.cbegin(), ak_sums.cend(), i)));
+  constexpr static auto A_k(std::uint16_t k, int_t r) -> float_t {
+    // boost forces floating point value return type, keep it throughout since
+    // for large r the number of points on faces dominate all others
+    const auto combinations = boost::math::binomial_coefficient<float_t>(d, k);
+    return combinations * std::pow(2, k) * std::pow(2 * r - 1, d - k);
   }
 
-  template <std::uniform_random_bit_generator RNG>
-  constexpr static auto random_sign(RNG &rng) -> int {
-    return 2 * std::uniform_int_distribution{0, 1}(rng)-1;
-  }
+
 };
 
 } // namespace lerw
diff --git a/include/lerw.hpp b/include/lerw.hpp
@@ -40,6 +40,10 @@ using PointType = typename PointTypeSelector<dim>::type;
 
 template <point P, Norm n> struct DirectionSelector;
 
+template <point P> struct DirectionSelector<P, Norm::L1> {
+  using type = L2Direction<P>;
+};
+
 template <point P> struct DirectionSelector<P, Norm::L2> {
   using type = L2Direction<P>;
 };
@@ -54,6 +58,10 @@ using DirectionType = typename DirectionSelector<P, norm>::type;
 // Length depends on Point because of the type used to store coordinates
 template <point P, Norm n> struct LengthSelector;
 
+template <point P> struct LengthSelector<P, Norm::L1> {
+  using type = Zipf<typename field<P>::type>;
+};
+
 template <point P> struct LengthSelector<P, Norm::L2> {
   using type = Pareto;
 };

diff --git a/plot.py b/plot.py
@@ -131,7 +131,7 @@ def plot_alpha_dependence(
 
     plt.xlabel("α")
     plt.ylabel("D")
-    plt.title(f"Fractal Dimension vs Alpha in {dimension}D LERW")
+    plt.title(f"Fractal Dimension vs Alpha LERW")
     plt.grid(True, alpha=0.3)
     plt.legend(loc="lower right")
 
@@ -205,8 +205,8 @@ def main():
     plt.figure(figsize=(10, 6))
     x_range = np.array([min(alpha_values), max(alpha_values)])
     plt.plot(x_range, x_range, "--", color="gray", alpha=0.7, label="α = D")
-    for dim, norm in product([2, 3], [Norm.L2, Norm.LINF]):
-        show_alpha_dependence(dim, norm, R_values, alpha_values, 1000)
+    for dim, norm in product([2, 3], [Norm.L1, Norm.L2, Norm.LINF]):
+        show_alpha_dependence(dim, norm, R_values, alpha_values, 500)
     plt.show()
 
 

diff --git a/src/main.cpp b/src/main.cpp
@@ -82,11 +82,18 @@ auto main(int argc, char *argv[]) -> int {
   auto computer = LERWComputer{[&seed_rng] { return std::mt19937{seed_rng()}; },
                                N, alpha, distance};
 
-  std::println(*out, "# D={}, R={}, N={}, α={}, Norm={}, seed={}", dimension,
-               distance, N, alpha, norm_to_string(norm), seed);
-
   const auto lengths = [&] {
     switch (switch_pair(dimension, norm)) {
+    case switch_pair(1, Norm::L1):
+      return computer.compute<1, Norm::L1>();
+    case switch_pair(2, Norm::L1):
+      return computer.compute<2, Norm::L1>();
+    case switch_pair(3, Norm::L1):
+      return computer.compute<3, Norm::L1>();
+    case switch_pair(4, Norm::L1):
+      return computer.compute<4, Norm::L1>();
+    case switch_pair(5, Norm::L1):
+      return computer.compute<5, Norm::L1>();
     case switch_pair(1, Norm::L2):
       return computer.compute<1, Norm::L2>();
     case switch_pair(2, Norm::L2):
@@ -97,8 +104,9 @@ auto main(int argc, char *argv[]) -> int {
       return computer.compute<4, Norm::L2>();
     case switch_pair(5, Norm::L2):
       return computer.compute<5, Norm::L2>();
-    case switch_pair(1, Norm::LINF):
-      return computer.compute<1, Norm::LINF>();
+      // TODO: LINF with d=1 is broken
+    // case switch_pair(1, Norm::LINF):
+    //   return computer.compute<1, Norm::LINF>();
     case switch_pair(2, Norm::LINF):
       return computer.compute<2, Norm::LINF>();
     case switch_pair(3, Norm::LINF):
@@ -112,6 +120,9 @@ auto main(int argc, char *argv[]) -> int {
     }
   }();
 
+  std::println(*out, "# D={}, R={}, N={}, α={}, Norm={}, seed={}", dimension,
+               distance, N, alpha, norm_to_string(norm), seed);
+
   for (auto l : lengths) {
     std::println(*out, "{}", l);
   }

diff --git a/tests/directions.cpp b/tests/directions.cpp
@@ -13,13 +13,10 @@
 
 using namespace lerw;
 
-using D2 = LinfDirection<ArrayPoint<2>>;
-using D3 = LinfDirection<ArrayPoint<3>>;
 using int_t = field<ArrayPoint<2>>::type;
 
-template <typename T>
-auto is_approximately_uniform(const std::vector<T> &vec,
-                              double tolerance) -> bool {
+template <class T>
+void check_approximately_uniform(const std::vector<T> &vec, double tolerance) {
   std::unordered_map<T, std::size_t> frequency;
   for (const T &element : vec) {
     frequency[element]++;
@@ -28,39 +25,65 @@ auto is_approximately_uniform(const std::vector<T> &vec,
   double expected_frequency =
       static_cast<double>(vec.size()) / frequency.size();
 
-  double max_allowed_deviation = expected_frequency * tolerance;
+  int max_allowed_deviation = expected_frequency * tolerance;
 
   for (const auto &pair : frequency) {
-    double deviation = std::abs(pair.second - expected_frequency);
-    if (deviation > max_allowed_deviation) {
-      return false;
-    }
-  }
+    int deviation = std::abs(pair.second - expected_frequency);
 
-  return true;
+    CHECK(deviation < max_allowed_deviation);
+  }
 }
 
-template <std::size_t Dim> void require_uniform(int_t r, std::size_t N) {
-  using P = ArrayPoint<Dim>;
+template <direction Direction>
+void check_uniform(Direction d, int_t r, std::size_t N) {
   auto rng = std::mt19937{};
-  auto points = std::vector<P>{N};
+  auto points = std::vector<typename Direction::result_type>{N};
   std::transform(points.cbegin(), points.cend(), points.begin(),
-                 [r, &rng](auto) { return LinfDirection<P>{}(r, rng); });
-  REQUIRE(is_approximately_uniform(points, 0.1));
+                 [r, &rng, &d](auto) { return d(r, rng); });
+  check_approximately_uniform(points, 0.1);
 }
 
-template <std::size_t Dim> void check_distance(int_t r, std::size_t N) {
-  using P = ArrayPoint<Dim>;
+template <direction Direction, Norm n>
+void require_length(int_t r, std::size_t N) {
   auto rng = std::mt19937{};
-  auto points = std::vector<P>{N};
+  auto points = std::vector<typename Direction::result_type>{N};
 
   std::ranges::for_each(
       std::ranges::iota_view(0) | std::views::take(N),
-      [r](auto p) { REQUIRE(norm<Norm::LINF>(p) == r); },
-      [r, &rng](auto) { return LinfDirection<ArrayPoint<Dim>>{}(r, rng); });
+      [r](auto p) { REQUIRE(norm<n>(p) == r); },
+      [r, &rng](auto) { return Direction{}(r, rng); });
+}
+
+TEST_CASE("L1") {
+  using D1 = L1Direction<ArrayPoint<1>>;
+  using D2 = L1Direction<ArrayPoint<2>>;
+  using D3 = L1Direction<ArrayPoint<3>>;
+  const std::size_t N = 100'000;
+
+  SECTION("uniform") {
+    check_uniform(D1{}, 1, N);
+    check_uniform(D1{}, 100, N);
+    check_uniform(D2{}, 1, N);
+    check_uniform(D2{}, 5, N);
+    check_uniform(D3{}, 1, N);
+    check_uniform(D3{}, 2, N);
+    check_uniform(L1Direction<ArrayPoint<4>>{}, 1, N);
+  }
+
+  SECTION("max r") {
+    const int_t r =
+        GENERATE(1, 2, 5, 100); // , std::numeric_limits<int_t>::max());
+    require_length<D1, Norm::L1>(r, N);
+    require_length<D2, Norm::L1>(r, N);
+    require_length<D3, Norm::L1>(r, N);
+    require_length<L1Direction<ArrayPoint<4>>, Norm::L1>(r, N);
+  }
 }
 
+// TODO: Broken for d=1, why?
 TEST_CASE("LINF") {
+  using D2 = LinfDirection<ArrayPoint<2>>;
+  using D3 = LinfDirection<ArrayPoint<3>>;
   SECTION("A_k") {
     // A square of sidelength 4 has 3 points per side
     // and 1 point per vertex
@@ -79,18 +102,22 @@ TEST_CASE("LINF") {
   }
 
   const std::size_t N = 100'000;
+
   SECTION("uniform") {
-    require_uniform<2>(1, N);
-    require_uniform<2>(5, N);
-    require_uniform<3>(1, N);
-    require_uniform<3>(2, N);
-    require_uniform<4>(1, N);
+    // check_uniform(LinfDirection<ArrayPoint<1>>{}, 1, N);
+    // require_uniform(LinfDirection<ArrayPoint<1>>{}, 100, N);
+    check_uniform(D2{}, 1, N);
+    check_uniform(D2{}, 5, N);
+    check_uniform(D3{}, 1, N);
+    check_uniform(D3{}, 2, N);
+    check_uniform(LinfDirection<ArrayPoint<4>>{}, 1, N);
   }
 
   SECTION("max r") {
     const int_t r = GENERATE(1, 2, 5, 100, std::numeric_limits<int_t>::max());
-    check_distance<2>(r, N);
-    check_distance<3>(r, N);
-    check_distance<4>(r, N);
+    // require_length<LinfDirection<ArrayPoint<1>>, Norm::LINF>(r, N);
+    require_length<D2, Norm::LINF>(r, N);
+    require_length<D3, Norm::LINF>(r, N);
+    require_length<LinfDirection<ArrayPoint<4>>, Norm::LINF>(r, N);
   }
 }