Print complete walks

CMakeLists.txt

@@ -11,10 +11,16 @@ include_directories(${CMAKE_SOURCE_DIR}/include)
 # including gtl as a system dependency prevents warnings when compiling it
 include_directories(SYSTEM ${CMAKE_SOURCE_DIR}/gtl/include)
 
-add_executable(lerw src/main.cpp)
+set(COMMON_SOURCES
+    src/config.cpp
+)
+
+add_executable(lerw_length src/main_length.cpp ${COMMON_SOURCES})
+add_executable(lerw_points src/main_points.cpp ${COMMON_SOURCES})
 
+foreach(target lerw_length lerw_points)
 # options from https://github.com/cpp-best-practices/cmake_template
-target_compile_options(lerw PUBLIC
+target_compile_options(${target} PUBLIC
 			    -Wall
   			    -Wextra # reasonable and standard
        			    -Wshadow # warn the user if a variable declaration shadows one from a parent context
@@ -41,20 +47,21 @@ target_compile_options(lerw PUBLIC
         		    -Wsuggest-override # warn if an overridden member function is not marked 'override' or 'final'
 			    )
 
-target_compile_options(lerw PUBLIC -O3 -march=native)
-target_compile_options(lerw PUBLIC -fconcepts-diagnostics-depth=4)
+target_compile_options(${target} PUBLIC -O3 -march=native)
+target_compile_options(${target} PUBLIC -fconcepts-diagnostics-depth=4)
 
 # silence warnings about portability of the gtl
 # https://gcc.gnu.org/onlinedocs/gcc-14.1.0/gcc/Warning-Options.html#index-Winterference-size
-target_compile_options(lerw PUBLIC -Wno-interference-size)
+target_compile_options(${target} PUBLIC -Wno-interference-size)
 
 # silence notes from the pstl-implementation
 # https://stackoverflow.com/a/23995391
-target_compile_options(lerw PUBLIC -fcompare-debug-second)
+target_compile_options(${target} PUBLIC -fcompare-debug-second)
 
-target_link_libraries(lerw PRIVATE
+target_link_libraries(${target} PRIVATE
   tbb
   boost_program_options)
+endforeach()
 
 # testing
 find_package(Catch2 3 REQUIRED)
@@ -76,4 +83,4 @@ include(Catch)
 catch_discover_tests(tests)
 
 # nix-build wants an 'install' target
-install(TARGETS lerw DESTINATION .)
+install(TARGETS lerw_length lerw_points DESTINATION .)

Makefile

@@ -10,7 +10,8 @@ build:
 
 install: build interface.py
 	mkdir -p $(INSTALL_DIR)/bin
-	cp result/lerw $(INSTALL_DIR)/bin
+	cp result/lerw_length $(INSTALL_DIR)/bin
+	cp result/lerw_points $(INSTALL_DIR)/bin
 	cp -n interface.py $(INSTALL_DIR)
 
 build_manual:

include/config.hpp

@@ -0,0 +1,23 @@
+#pragma once
+
+#include <cstddef>
+#include <string>
+
+#include "utils.hpp"
+
+namespace lerw {
+
+struct Configuration {
+  Norm norm;
+  std::size_t dimension;
+  std::size_t num_samples;
+  double distance;
+  double alpha;
+  std::string output_path;
+  std::size_t seed;
+  bool help_requested;
+};
+
+Configuration parse_command_line(int argc, char *argv[]);
+
+} // namespace lerw

include/generator.hpp

@@ -8,8 +8,8 @@
 namespace lerw {
 
 template <stopper Stopper, stepper Stepper> struct RandomWalkGenerator {
-  Stopper stopper;
   Stepper stepper;
+  Stopper stopper;
 
   template <std::uniform_random_bit_generator RNG>
   constexpr auto operator()(RNG &rng) -> auto {
@@ -24,8 +24,8 @@ template <stopper Stopper, stepper Stepper> struct RandomWalkGenerator {
 
 template <stopper Stopper, stepper Stepper>
 struct LoopErasedRandomWalkGenerator {
-  Stopper stopper;
   Stepper stepper;
+  Stopper stopper;
 
   template <std::uniform_random_bit_generator RNG>
   constexpr auto operator()(RNG &rng) -> auto {

include/lerw.hpp

@@ -73,26 +73,43 @@ template <point P> struct LengthSelector<P, Norm::LINF> {
 template <point P, Norm n>
 using LengthType = typename LengthSelector<P, n>::type;
 
-struct LERWComputer {
+struct DistanceLerwComputer {
   std::function<std::mt19937()> rng_factory;
   std::size_t N;
   double alpha;
   double distance;
-  template <std::size_t dim, Norm norm> auto compute() const {
+  template <std::size_t dim, Norm norm, class Projection = std::identity>
+  auto compute(Projection projection = {}) const {
     using point_t = PointType<dim>;
-    return compute_lerw_lengths(
-        [alpha = alpha]() {
-          return LDStepper{LengthType<point_t, norm>{alpha}, DirectionType<point_t, norm>{}};
+    return compute_lengths(
+        [alpha = alpha, distance = distance] {
+          return LoopErasedRandomWalkGenerator{
+              LDStepper{LengthType<point_t, norm>{alpha},
+                        DirectionType<point_t, norm>{}},
+              DistanceStopper<norm>{distance}};
         },
-        [distance = distance]() { return DistanceStopper<norm>{distance}; },
-        rng_factory, N);
+        rng_factory, projection, N);
   }
 };
 
-template <class GeneratorFactory, class RNGFactory>
+struct LengthLerwComputer {
+  double alpha;
+  size_t length;
+
+  template <std::size_t dim, Norm norm, std::uniform_random_bit_generator RNG>
+  auto compute(RNG &rng) const {
+    using point_t = PointType<dim>;
+    return LoopErasedRandomWalkGenerator{
+        LDStepper{LengthType<point_t, norm>{alpha},
+                  DirectionType<point_t, norm>{}},
+        LengthStopper{length}}(rng);
+  }
+};
+
+template <class GeneratorFactory, class RNGFactory, class Projection>
 auto compute_lengths(GeneratorFactory &&generator_factory,
-                     RNGFactory &&rng_factory,
-                     size_t N) -> std::vector<size_t> {
+                     RNGFactory &&rng_factory, Projection projection,
+                     size_t N) {
   auto generators = std::vector<decltype(generator_factory())>{};
   auto rngs = std::vector<decltype(rng_factory())>{};
 
@@ -105,52 +122,13 @@ auto compute_lengths(GeneratorFactory &&generator_factory,
 
   std::vector<size_t> lengths(N);
 
-  std::transform(
-      std::execution::par_unseq, generators.begin(), generators.end(),
-      rngs.begin(), lengths.begin(),
-      [](auto generator, auto rng) { return generator(rng).size(); });
+  std::transform(std::execution::par_unseq, generators.begin(),
+                 generators.end(), rngs.begin(), lengths.begin(),
+                 [projection](auto generator, auto rng) {
+                   return projection(generator(rng));
+                 });
 
   return lengths;
 }
 
-template <class StepperFactory, class StopperFactory, class RNGFactory>
-auto compute_lerw_lengths(StepperFactory &&stepper_factory,
-                          StopperFactory &&stopper_factory,
-                          RNGFactory &&rng_factory,
-                          std::size_t n_samples) -> auto {
-  auto generator_factory = [&stopper_factory, &stepper_factory] {
-    return LoopErasedRandomWalkGenerator{stopper_factory(), stepper_factory()};
-  };
-  return compute_lengths(generator_factory, rng_factory, n_samples);
-}
-
-template <class GeneratorFactory, class RNGFactory>
-auto compute_average_length(GeneratorFactory &&generator_factory,
-                            RNGFactory &&rng_factory, size_t N) -> double {
-  assert(N != 0);
-  return std::ranges::fold_left_first(
-             compute_lengths(std::move(generator_factory),
-                             std::move(rng_factory), N),
-             std::plus{}) /
-         static_cast<double>(N);
-}
-
-template <class StepperFactory, class RNGFactory>
-auto compute_lerw_average_lengths(StepperFactory &&stepper_factory,
-                                  RNGFactory &&rng_factory,
-                                  const std::vector<double> &distances,
-                                  std::size_t n_samples) -> auto {
-  auto results = std::vector<std::pair<double, double>>{};
-  for (const auto &d : distances) {
-    auto stopper_factory = [d] { return DistanceStopper<Norm::L2>{d}; };
-    auto l = compute_average_length(
-        [&stopper_factory, &stepper_factory] {
-          return LoopErasedRandomWalkGenerator{stopper_factory(),
-                                               stepper_factory()};
-        },
-        rng_factory, n_samples);
-    results.emplace_back(d, l);
-  }
-  return results;
-}
 } // namespace lerw

include/utils.hpp

@@ -80,4 +80,10 @@ template <Norm N, class... T> constexpr auto norm(T... args) -> auto {
   return norm_selector<N>{}(args...);
 }
 
+
+// helper for switch
+constexpr auto switch_pair(std::size_t dimension, Norm norm) -> std::size_t {
+  return (dimension << 2) + static_cast<size_t>(norm);
+}
+
 } // namespace lerw

interface.py

@@ -7,7 +7,8 @@
 
 # relative to the location where this is called from
 DATA_DIR: Path = Path("data")
-CPP_EXECUTABLE: Path = Path("bin") / "lerw"
+CPP_LENGTHS_EXE: Path = Path("bin") / "lerw_length"
+CPP_POINTS_EXE: Path = Path("bin") / "lerw_points"
 
 
 class Norm(Enum):
@@ -16,6 +17,71 @@ class Norm(Enum):
     L2 = 2
 
 
+def get_walk(
+    number_of_steps: int,
+    alpha: float,
+    norm: Norm,
+    seed: int = 3,
+    recompute: bool = False,
+) -> npt.NDArray[np.int64]:
+    """Get Points of a single walk.
+
+    Arguments match output of "<cpp_exe> --help".
+    """
+    DATA_DIR.mkdir(exist_ok=True)
+    dimension = 2
+    distance = 0.0  # ignored by cpp
+
+    filename = _format_filename(
+        "walk", dimension, distance, number_of_steps, alpha, norm, seed
+    )
+    file_path = DATA_DIR / filename
+
+    if recompute:
+        file_path.unlink(missing_ok=True)
+
+    if not file_path.exists():
+        cmd = [
+            Path.cwd() / CPP_POINTS_EXE,
+            "--dimension",
+            dimension,
+            "--distance",
+            distance,
+            "--number_of_walks",
+            number_of_steps,
+            "--alpha",
+            alpha,
+            "--norm",
+            norm.name,
+            "--output",
+            file_path,
+            "--seed",
+            seed,
+        ]
+
+        result = subprocess.run(
+            list(map(str, cmd)),
+            capture_output=True,
+            text=True,
+            check=False,
+        )
+
+        # Check for any output, which indicates an error
+        if result.stdout or result.stderr or result.returncode != 0:
+            print(
+                f"Failed to call C++:\n{result.stderr}\n\n{result.stdout}",
+                file=sys.stderr,
+            )
+            raise subprocess.CalledProcessError(
+                returncode=result.returncode or 1,
+                cmd=cmd,
+                output=result.stdout,
+                stderr=result.stderr,
+            )
+
+    return np.genfromtxt(file_path, dtype=np.int64, comments="#", delimiter=",")
+
+
 def get_walk_lengths(
     dimension: int,
     distance: float,
@@ -31,15 +97,17 @@ def get_walk_lengths(
     """
     DATA_DIR.mkdir(exist_ok=True)
 
-    filename = _format_filename(dimension, distance, number_of_walks, alpha, norm, seed)
+    filename = _format_filename(
+        "walks", dimension, distance, number_of_walks, alpha, norm, seed
+    )
     file_path = DATA_DIR / filename
 
     if recompute:
         file_path.unlink(missing_ok=True)
 
     if not file_path.exists():
         cmd = [
-            Path.cwd() / CPP_EXECUTABLE,
+            Path.cwd() / CPP_LENGTHS_EXE,
             "--dimension",
             dimension,
             "--distance",
@@ -80,17 +148,32 @@ def get_walk_lengths(
 
 
 def _format_filename(
+    prefix: str,
     dimension: int,
     distance: float,
     number_of_walks: int,
     alpha: float,
     norm: Norm,
     seed: int = 42,
 ) -> str:
-    return f"walks_dim{dimension}_dist{distance}_n{number_of_walks}_a{alpha}_{norm.name}_rng{seed}.txt"
+    return f"{prefix}_dim{dimension}_dist{distance}_n{number_of_walks}_a{alpha}_{norm.name}_rng{seed}.txt"
 
 
 def test():
+    num_steps = 10
+    args = {
+        "alpha": 0.5,
+        "norm": Norm.L2,
+        "seed": 2,
+    }
+    file = Path(DATA_DIR) / _format_filename("walk", 2, 0.0, num_steps, **args)
+    file.unlink(missing_ok=True)
+
+    walk = get_walk(num_steps, **args)
+    assert file.exists()
+    assert walk.shape == (11, 2)
+    assert all(walk[10, :] == [16, -1077])
+
     args = {
         "dimension": 2,
         "distance": 5000,
@@ -99,7 +182,7 @@ def test():
         "norm": Norm.L2,
         "seed": 2,
     }
-    file = Path(DATA_DIR) / _format_filename(**args)
+    file = Path(DATA_DIR) / _format_filename("walks", **args)
     file.unlink(missing_ok=True)
 
     walks = get_walk_lengths(**args)