Merge pull request #5 from mino2357/speed-up

高速化

.gitignore

@@ -16,4 +16,7 @@ Cargo.lock
 *.png
 *.log
 
+*.svg
+perf.*
+
 !images/*.png

src/gen_grid2d.rs

@@ -14,8 +14,8 @@ impl Grid2D {
 
     #[allow(dead_code)]
     pub fn distance_square(&self, point: &Grid2D) -> f64 {
-        let dx2 = (self.x - point.x).powf(2.0);
-        let dy2 = (self.y - point.y).powf(2.0);
+        let dx2 = (self.x - point.x) * (self.x - point.x);
+        let dy2 = (self.y - point.y) * (self.y - point.y);
         dx2 + dy2
     }
 }
@@ -104,38 +104,53 @@ impl Points2D {
         radius: f64,
     ) -> Grid2D {
         let epsilon: f64 = 1.0;
-        let sigma: f64 = 0.9 * (2.0_f64).powf(-1.0 / 6.0) / (self.points.len() as f64).powf(0.5);
+        let sigma: f64 = 0.9 * (2.0_f64).powf(-1.0 / 6.0) / (self.points.len() as f64).sqrt();
         let mut f_x = 0.0;
         let mut f_y = 0.0;
+        let sigma2 = sigma * sigma;
+        let sigma4 = sigma2 * sigma2;
+        let sigma8 = sigma4 * sigma4;
+        let sigma6 = sigma2 * sigma4;
+        let sigma12 = sigma8 * sigma4;
         let near = tree.neighbor_search(&self.points[index], radius);
         for k in near.iter() {
             if index != *k {
                 let dx = self.points[index].x - self.points[*k].x;
                 let dy = self.points[index].y - self.points[*k].y;
-                let r = (dx.powf(2.0) + dy.powf(2.0)).powf(0.5);
+                let r = (dx * dx + dy * dy).sqrt();
+                let r2 = r * r;
+                let r4 = r2 * r2;
+                let r8 = r4 * r4;
+                let r13 = r8 * r4 * r;
+                let r7 = r2 * r4 * r;
                 f_x += 4.0
                     * epsilon
-                    * (12.0 * sigma.powf(12.0) / r.powf(13.0)
-                        - 6.0 * sigma.powf(6.0) / r.powf(7.0))
+                    * (12.0 * sigma12 / r13
+                        - 6.0 * sigma6 / r7)
                     * dx;
                 f_y += 4.0
                     * epsilon
-                    * (12.0 * sigma.powf(12.0) / r.powf(13.0)
-                        - 6.0 * sigma.powf(6.0) / r.powf(7.0))
+                    * (12.0 * sigma12 / r13
+                        - 6.0 * sigma6 / r7)
                     * dy;
             }
         }
         for i in 0..boundary.points.len() {
             let dx = self.points[index].x - boundary.points[i].x;
             let dy = self.points[index].y - boundary.points[i].y;
-            let r = (dx.powf(2.0) + dy.powf(2.0)).powf(0.5);
+            let r = (dx * dx + dy * dy).sqrt();
+            let r2 = r * r;
+            let r4 = r2 * r2;
+            let r8 = r4 * r4;
+            let r13 = r8 * r4 * r;
+            let r7 = r2 * r4 * r;
             f_x += 4.0
                 * epsilon
-                * (12.0 * sigma.powf(12.0) / r.powf(13.0) - 6.0 * sigma.powf(6.0) / r.powf(7.0))
+                * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7)
                 * dx;
             f_y += 4.0
                 * epsilon
-                * (12.0 * sigma.powf(12.0) / r.powf(13.0) - 6.0 * sigma.powf(6.0) / r.powf(7.0))
+                * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7)
                 * dy;
         }
         Grid2D { x: f_x, y: f_y }

src/main.rs

@@ -12,7 +12,7 @@ mod kd_tree;
 fn draw_graph(i: usize, vec: &gen_grid2d::Points2D, boundary: &gen_grid2d::Points2D) {
     let out_file_name = format!("{:04}", i).to_string() + ".png";
 
-    let root = BitMapBackend::new(&out_file_name, (1080, 1080)).into_drawing_area();
+    let root = BitMapBackend::new(&out_file_name, (2000, 2000)).into_drawing_area();
 
     root.fill(&WHITE).unwrap();
 
@@ -29,7 +29,7 @@ fn draw_graph(i: usize, vec: &gen_grid2d::Points2D, boundary: &gen_grid2d::Point
     chart
         .draw_series(PointSeries::of_element(
             (0..vec.points.len()).map(|i| (vec.points[i].x, vec.points[i].y)),
-            2,
+            1,
             ShapeStyle::from(&RED).filled(),
             &|coord, size, style| EmptyElement::at(coord) + Circle::new((0, 0), size, style),
         ))
@@ -38,7 +38,7 @@ fn draw_graph(i: usize, vec: &gen_grid2d::Points2D, boundary: &gen_grid2d::Point
     chart
         .draw_series(PointSeries::of_element(
             (0..boundary.points.len()).map(|i| (boundary.points[i].x, boundary.points[i].y)),
-            2,
+            1,
             ShapeStyle::from(&BLUE).filled(),
             &|coord, size, style| EmptyElement::at(coord) + Circle::new((0, 0), size, style),
         ))
@@ -90,8 +90,8 @@ fn main() {
     let seed: [u8; 32] = [1; 32];
     let mut rng: rand::rngs::StdRng = rand::SeedableRng::from_seed(seed);
 
-    let num_point: usize = 400;
-    let num_boundary: usize = (std::f64::consts::PI * (num_point as f64).powf(0.5)) as usize;
+    let num_point: usize = 100000;
+    let num_boundary: usize = (std::f64::consts::PI * (num_point as f64).sqrt()) as usize;
 
     let mut vec = gen_grid2d::Points2D::new();
     let mut boundary = gen_grid2d::Points2D::new();
@@ -107,10 +107,10 @@ fn main() {
         }
     }
 
-    let radius = 4.0 * std::f64::consts::PI / (num_boundary as f64);
+    let radius = 6.0 * std::f64::consts::PI / (num_boundary as f64);
     println!("radius = {}", radius);
 
-    let max_counter = 50;
+    let max_counter = 500;
 
     let mut tree = kd_tree::KDTree::construct_kd_tree(&vec);
 
@@ -121,7 +121,7 @@ fn main() {
 
     for i in 0..max_counter {
         draw_graph(i, &vec, &boundary);
-        for _ in 0..100 {
+        for _ in 0..10 {
             vec.euler_step_by_near_points(&boundary, &tree, radius);
             //vec.euler_step(&boundary);
             tree = kd_tree::KDTree::construct_kd_tree(&vec);