Done

.gitignore

@@ -1 +1,4 @@
 /target
+*.csv
+*.pdf
+*.png

Cargo.toml

@@ -1,14 +1,24 @@
+[workspace]
+members = [
+  "crates/*",
+]
+
+[workspace.package]
+version = "0.0.1"
+authors = ["Mo Bitar <mo8it@proton.me>"]
+edition = "2021"
+license = "AGPL-3.0"
+
 [package]
 name = "pivot-saw"
 description = "Self-avoiding walk with the pivot algorithm"
-version = "0.1.0"
-edition = "2021"
-license = "AGPL-3.0"
-repository = "https://codeberg.org/mo8it/pivot-saw"
+version.workspace = true
+authors.workspace = true
+edition.workspace = true
+license.workspace = true
 
 [dependencies]
-bevy = "0.10.1"
-nalgebra = { version = "0.32.2", default-features = false, features = ["std"] }
+nalgebra = { version = "0.32.3", default-features = false, features = ["std"] }
 rand = { version = "0.8.5", features = ["small_rng"] }
 
 [profile.dev]

crates/analysis/Cargo.toml

@@ -0,0 +1,11 @@
+[package]
+name = "analysis"
+publish = false
+version.workspace = true
+authors.workspace = true
+edition.workspace = true
+license.workspace = true
+
+[dependencies]
+pivot-saw = { path = "../.." }
+rayon = "1.7.0"

crates/analysis/Project.toml

@@ -0,0 +1,4 @@
+[deps]
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

crates/analysis/main.jl

@@ -0,0 +1,30 @@
+using DelimitedFiles
+using Plots
+using LsqFit
+
+@. model(N, p) = N^(2 * p[1])
+
+function main()
+    data = readdlm("analysis.csv", ',')
+
+    ns = view(data, :, 1)
+    squared_r_ees = view(data, :, 2)
+
+    expected_ν = 0.588
+    expected(N) = N^(2 * expected_ν)
+
+    p0 = [expected_ν]
+    fit = curve_fit(model, ns, squared_r_ees, p0)
+    simulation_ν = fit.param[1]
+
+    p = plot(; xlabel="N", ylabel="<Rₑₑ²>", xscale=:log10, yscale=:log10, title="10⁵ steps", dpi=300)
+
+    plot!(p, ns, squared_r_ees; label="Simulation ν = $(round(simulation_ν; digits=3))")
+    plot!(p, ns, map(N -> expected(N), ns); label="Expected ν = $expected_ν")
+
+    savefig(p, "analysis.png")
+
+    return nothing
+end
+
+main()

crates/analysis/src/main.rs

@@ -0,0 +1,49 @@
+use rayon::prelude::*;
+use std::{
+    fs::OpenOptions,
+    io::{self, Write},
+};
+
+use pivot_saw::Walk;
+
+fn main() -> Result<(), io::Error> {
+    let n_range = 100..6_001;
+    let step_by = 200;
+    let max_n_steps = 100_000;
+
+    let mut file = OpenOptions::new()
+        .create(true)
+        .write(true)
+        .truncate(true)
+        .open("analysis.csv")?;
+
+    let squared_r_ees = n_range
+        .clone()
+        .into_par_iter()
+        .step_by(step_by)
+        .map(|n| {
+            let mut n_steps = 1_u64;
+            let mut mean_r_ee_squared = (n * n) as f64;
+
+            let mut walk = Walk::new(n);
+
+            for _ in 0..max_n_steps {
+                walk.pivot();
+
+                n_steps += 1;
+                let r_ee = walk.points.last().unwrap();
+                let r_ee_squared = r_ee.x * r_ee.x + r_ee.y * r_ee.y;
+                mean_r_ee_squared += (r_ee_squared as f64 - mean_r_ee_squared) / n_steps as f64;
+            }
+
+            mean_r_ee_squared
+        })
+        .collect::<Vec<_>>();
+
+    for (n, squared_r_ee) in n_range.step_by(step_by).zip(squared_r_ees) {
+        writeln!(file, "{n},{squared_r_ee}")?;
+    }
+    file.flush()?;
+
+    Ok(())
+}

crates/visualization/Cargo.toml

@@ -0,0 +1,11 @@
+[package]
+name = "visualization"
+publish = false
+version.workspace = true
+authors.workspace = true
+edition.workspace = true
+license.workspace = true
+
+[dependencies]
+bevy = "0.11.0"
+pivot-saw = { path = "../.." }

crates/visualization/src/main.rs

@@ -0,0 +1,60 @@
+use bevy::{
+    prelude::{
+        App, Assets, Camera3dBundle, Color, Commands, DefaultPlugins, Mesh, PbrBundle, ResMut,
+        Resource, StandardMaterial, Startup, Transform, Update, Vec3,
+    },
+    render::render_resource::PrimitiveTopology,
+};
+
+use pivot_saw::Walk;
+
+#[derive(Resource)]
+struct WalkRes(Walk);
+
+fn setup(mut commands: Commands) {
+    commands.spawn(Camera3dBundle {
+        transform: Transform::from_xyz(80.0, 80.0, 80.0)
+            .looking_at(Vec3::new(0.0, 0.0, 0.0), Vec3::Y),
+        ..Default::default()
+    });
+}
+
+fn pivot(
+    mut commands: Commands,
+    mut meshes: ResMut<Assets<Mesh>>,
+    mut materials: ResMut<Assets<StandardMaterial>>,
+    mut walk: ResMut<WalkRes>,
+) {
+    let blue = materials.add(Color::BLUE.into());
+
+    let mut mesh = Mesh::new(PrimitiveTopology::LineStrip);
+    mesh.insert_attribute(
+        Mesh::ATTRIBUTE_POSITION,
+        walk.0
+            .points
+            .iter()
+            .map(|p| Vec3::new(p.x as f32, p.y as f32, p.z as f32))
+            .collect::<Vec<_>>(),
+    );
+
+    commands.spawn(PbrBundle {
+        mesh: meshes.add(mesh),
+        material: blue,
+        ..Default::default()
+    });
+
+    walk.0.pivot();
+}
+
+fn main() {
+    let n = 500;
+
+    let walk = WalkRes(Walk::new(n));
+
+    App::new()
+        .add_plugins(DefaultPlugins)
+        .insert_resource(walk)
+        .add_systems(Startup, setup)
+        .add_systems(Update, pivot)
+        .run();
+}

src/lib.rs

@@ -1,10 +1,3 @@
-use bevy::{
-    prelude::{
-        App, Assets, Camera3dBundle, Color, Commands, DefaultPlugins, Mesh, PbrBundle, ResMut,
-        Resource, StandardMaterial, Transform, Vec3,
-    },
-    render::render_resource::PrimitiveTopology,
-};
 use nalgebra::base::{Matrix3, Vector3};
 use rand::{distributions::Uniform, prelude::Distribution, rngs::SmallRng, Rng, SeedableRng};
 use std::collections::HashSet;
@@ -17,8 +10,8 @@ struct Transformations {
     dist: Uniform<usize>,
 }
 
-impl Transformations {
-    pub fn new() -> Self {
+impl Default for Transformations {
+    fn default() -> Self {
         let len = 23;
         let mut ms = Vec::with_capacity(len);
 
@@ -28,40 +21,41 @@ impl Transformations {
             for v1 in [1, -1] {
                 for v2 in [1, -1] {
                     for ind in 0..3 {
-                        if unit_matrix_skipped {
-                            let mut m = Matrix3::zeros();
+                        for shift in 0..2 {
+                            if unit_matrix_skipped {
+                                let mut m = Matrix3::zeros();
 
-                            *m.index_mut((0, ind)) = v0;
-                            *m.index_mut((1, (ind + 1) % 3)) = v1;
-                            *m.index_mut((2, (ind + 2) % 3)) = v2;
+                                *m.index_mut((0, ind)) = v0;
+                                *m.index_mut((1, (ind + shift + 1) % 3)) = v1;
+                                *m.index_mut((2, (ind + shift + 2) % 3)) = v2;
 
-                            ms.push(m);
-                        } else {
-                            unit_matrix_skipped = true;
+                                ms.push(m);
+                            } else {
+                                unit_matrix_skipped = true;
+                            }
                         }
                     }
                 }
             }
         }
 
-        debug_assert_eq!(len, ms.len());
-
         Self {
             ms,
             dist: Uniform::new(0, len),
         }
     }
+}
 
-    pub fn rand(&self, rng: &mut impl Rng) -> &Transformation {
+impl Transformations {
+    fn rand(&self, rng: &mut impl Rng) -> &Transformation {
         let ind = self.dist.sample(rng);
 
         unsafe { self.ms.get_unchecked(ind) }
     }
 }
 
-#[derive(Resource)]
-struct Walk {
-    points: Vec<Point>,
+pub struct Walk {
+    pub points: Vec<Point>,
     new_points: Vec<Point>,
     pivot_dist: Uniform<usize>,
     set: HashSet<Point>,
@@ -83,7 +77,7 @@ impl Walk {
 
         let rng = SmallRng::seed_from_u64(42);
 
-        let transformations = Transformations::new();
+        let transformations = Transformations::default();
 
         Self {
             points,
@@ -108,7 +102,7 @@ impl Walk {
             let pivot = unsafe { s1.get_unchecked(pivot_ind) };
 
             for point in s2.iter() {
-                self.new_points.push((m * (point - pivot)) + pivot);
+                self.new_points.push(m * (point - pivot) + pivot);
             }
 
             if s1
@@ -125,50 +119,3 @@ impl Walk {
         }
     }
 }
-
-fn setup(mut commands: Commands) {
-    commands.spawn(Camera3dBundle {
-        transform: Transform::from_xyz(80.0, 80.0, 80.0)
-            .looking_at(Vec3::new(0.0, 0.0, 0.0), Vec3::Y),
-        ..Default::default()
-    });
-}
-
-fn pivot(
-    mut commands: Commands,
-    mut meshes: ResMut<Assets<Mesh>>,
-    mut materials: ResMut<Assets<StandardMaterial>>,
-    mut walk: ResMut<Walk>,
-) {
-    let blue = materials.add(Color::BLUE.into());
-
-    let mut mesh = Mesh::new(PrimitiveTopology::LineStrip);
-    mesh.insert_attribute(
-        Mesh::ATTRIBUTE_POSITION,
-        walk.points
-            .iter()
-            .map(|p| Vec3::new(p.x as f32, p.y as f32, p.z as f32))
-            .collect::<Vec<_>>(),
-    );
-
-    commands.spawn(PbrBundle {
-        mesh: meshes.add(mesh),
-        material: blue,
-        ..Default::default()
-    });
-
-    walk.pivot();
-}
-
-fn main() {
-    let n = 250;
-
-    let walk = Walk::new(n);
-
-    App::new()
-        .add_plugins(DefaultPlugins)
-        .insert_resource(walk)
-        .add_startup_system(setup)
-        .add_system(pivot)
-        .run();
-}