diff --git a/Makefile b/Makefile index a84f2a3..9ec615b 100644 --- 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 index d49f32a..6b7a75c 100644 --- a/include/directions.hpp +++ b/include/directions.hpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -31,6 +32,95 @@ template struct L2Direction { } }; +template +constexpr auto random_sign(RNG &rng) -> std::int8_t { + using i = std::int8_t; + auto d = std::uniform_int_distribution{0, 1}; + return static_cast(2 * d(rng) - 1); +} + +// uniformly pick direction from surface of L1-ball +template 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::type; + + static constexpr std::size_t d = dim(); + static_assert(d > 0, "Out-of-bounds stuff happens when d = 0."); + + template + 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{}(proposed.cbegin(), proposed.cend()); + } + + template + constexpr static auto get_integer_solution(int_t r, + RNG &rng) -> std::array { + auto bars = sample_iota(r, rng); + bars[d - 1] = r + d; + std::sort(bars.begin(), bars.end()); + auto stars = std::array{}; + 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 + constexpr static auto reject_for_zeros(InputIt first, InputIt last, + RNG &rng) -> bool { + auto d = std::uniform_int_distribution{0, 1}; + return std::transform_reduce( + first, last, false, std::logical_or{}, + [&d, &rng](auto c) { return c == 0 && d(rng) == 0; }); + } + + template + constexpr static auto sample_iota(int_t r, RNG &rng) -> std::array { + // equivalent to + // auto indices = std::vector(r + d - 1); + // std::iota(indices.begin(), indices.end(), 1); + // auto result = std::array{}; + // 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{1, r + d - 1}; + auto result = std::array{}; + 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 struct LinfDirection { @@ -53,26 +143,28 @@ template struct LinfDirection { template constexpr auto operator()(int_t r, RNG &rng) -> Point { - r = std::abs(r); const auto k = choose_k(r, rng); auto coordinates = std::array{}; 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{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{}(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(d, k); - return combinations * std::pow(2, k) * std::pow(2 * r - 1, d - k); + template + 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{0, ak_sums[d - 1]}(rng); + return static_cast( + 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 { @@ -88,20 +180,14 @@ template struct LinfDirection { return ak_sums; } - template - 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{0, ak_sums[d - 1]}(rng); - return static_cast( - 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(d, k); + return combinations * std::pow(2, k) * std::pow(2 * r - 1, d - k); } - template - 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 index d99c134..6e7a943 100644 --- a/include/lerw.hpp +++ b/include/lerw.hpp @@ -40,6 +40,10 @@ using PointType = typename PointTypeSelector::type; template struct DirectionSelector; +template struct DirectionSelector { + using type = L2Direction

; +}; + template struct DirectionSelector { using type = L2Direction

; }; @@ -54,6 +58,10 @@ using DirectionType = typename DirectionSelector::type; // Length depends on Point because of the type used to store coordinates template struct LengthSelector; +template struct LengthSelector { + using type = Zipf::type>; +}; + template struct LengthSelector { using type = Pareto; }; diff --git a/plot.py b/plot.py index 7a3c9d3..2ca06e7 100644 --- 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 index 61ee35b..7d94561 100644 --- 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 index 8a48908..88bbfc1 100644 --- a/tests/directions.cpp +++ b/tests/directions.cpp @@ -13,13 +13,10 @@ using namespace lerw; -using D2 = LinfDirection>; -using D3 = LinfDirection>; using int_t = field>::type; -template -auto is_approximately_uniform(const std::vector &vec, - double tolerance) -> bool { +template +void check_approximately_uniform(const std::vector &vec, double tolerance) { std::unordered_map frequency; for (const T &element : vec) { frequency[element]++; @@ -28,39 +25,65 @@ auto is_approximately_uniform(const std::vector &vec, double expected_frequency = static_cast(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 void require_uniform(int_t r, std::size_t N) { - using P = ArrayPoint; +template +void check_uniform(Direction d, int_t r, std::size_t N) { auto rng = std::mt19937{}; - auto points = std::vector

{N}; + auto points = std::vector{N}; std::transform(points.cbegin(), points.cend(), points.begin(), - [r, &rng](auto) { return LinfDirection

{}(r, rng); }); - REQUIRE(is_approximately_uniform(points, 0.1)); + [r, &rng, &d](auto) { return d(r, rng); }); + check_approximately_uniform(points, 0.1); } -template void check_distance(int_t r, std::size_t N) { - using P = ArrayPoint; +template +void require_length(int_t r, std::size_t N) { auto rng = std::mt19937{}; - auto points = std::vector

{N}; + auto points = std::vector{N}; std::ranges::for_each( std::ranges::iota_view(0) | std::views::take(N), - [r](auto p) { REQUIRE(norm(p) == r); }, - [r, &rng](auto) { return LinfDirection>{}(r, rng); }); + [r](auto p) { REQUIRE(norm(p) == r); }, + [r, &rng](auto) { return Direction{}(r, rng); }); +} + +TEST_CASE("L1") { + using D1 = L1Direction>; + using D2 = L1Direction>; + using D3 = L1Direction>; + 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>{}, 1, N); + } + + SECTION("max r") { + const int_t r = + GENERATE(1, 2, 5, 100); // , std::numeric_limits::max()); + require_length(r, N); + require_length(r, N); + require_length(r, N); + require_length>, Norm::L1>(r, N); + } } +// TODO: Broken for d=1, why? TEST_CASE("LINF") { + using D2 = LinfDirection>; + using D3 = LinfDirection>; 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>{}, 1, N); + // require_uniform(LinfDirection>{}, 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>{}, 1, N); } SECTION("max r") { const int_t r = GENERATE(1, 2, 5, 100, std::numeric_limits::max()); - check_distance<2>(r, N); - check_distance<3>(r, N); - check_distance<4>(r, N); + // require_length>, Norm::LINF>(r, N); + require_length(r, N); + require_length(r, N); + require_length>, Norm::LINF>(r, N); } }