diff --git a/Cargo.toml b/Cargo.toml index a07da8c..4d7059e 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,5 +1,5 @@ [package] -name = "MonteCarlo-Rust" +name = "montecarlo-rust" version = "0.1.0" edition = "2021" diff --git a/src/bounds.rs b/src/bounds.rs index b375edb..53f42cb 100644 --- a/src/bounds.rs +++ b/src/bounds.rs @@ -12,19 +12,15 @@ impl Bounds { return Bounds {x, y}; } - pub fn get_area(self) -> f64 { + pub fn get_area(&self) -> f64 { return self.x.get_length() * self.y.get_length(); } - pub fn get_random_point(self) -> Point { + pub fn get_random_point(&self) -> Point { return Point::new(self.x.get_random_value(), self.y.get_random_value()); } - pub fn get_x(self) -> LinearBounds { + pub fn get_x(&self) -> LinearBounds { return self.x; } - - pub fn get_y(self) -> LinearBounds { - return self.y; - } } diff --git a/src/function.rs b/src/function.rs new file mode 100644 index 0000000..35f24f5 --- /dev/null +++ b/src/function.rs @@ -0,0 +1,61 @@ +use std::sync::Arc; +use std::thread::JoinHandle; +use crate::bounds::Bounds; +use crate::linear_bounds::LinearBounds; +use crate::point::Point; + +#[derive(Clone, Copy)] +pub struct Function{} + +fn f(x: &f64) -> f64 { + return x.powi(4) - 4_f64 * x.powi(3) + 18_f64 * x.powi(2) - 12_f64 * x - 69_f64 +} + +fn big_f(x: &f64) -> f64 { + return 0.2 * x.powi(5) - x.powi(4) + 6_f64 * x.powi(3) - 6_f64 * x.powi(2) - 69_f64 * x; +} + +impl Function { + pub fn new() -> Function { + return Function {}; + } + + pub fn integrate(&self, bounds: &LinearBounds) -> f64 { + return big_f(&bounds.get_higher()) - big_f(&bounds.get_lower()); + } + + fn includes(&self, point: &Point) -> bool { + let y_0: f64 = f(&point.get_x()); + return y_0.abs() > point.get_y().abs() && y_0 * point.get_y() >= 0_f64; + } + + pub fn approximate(&self, bounds: &Bounds, samples: usize, thread_cnt: usize) -> f64 { + let mut hits: u64 = 0; + let mut threads: Vec> = vec![]; + + let function: Arc = Arc::new(*self); + let bounds: Arc = Arc::new(*bounds); + + for _ in 0..thread_cnt { + let function = function.clone(); + let bounds = bounds.clone(); + + threads.push(std::thread::spawn(move || { + let mut hits: u64 = 0; + for _ in 0..samples { + let point: Point = bounds.get_random_point(); + if function.includes(&point) { + hits += 1; + } + } + return hits; + })) + } + + for thread in threads { + hits += thread.join().unwrap(); + } + + return (hits as f64 / (thread_cnt * samples) as f64) * bounds.get_area(); + } +} \ No newline at end of file diff --git a/src/linear_bounds.rs b/src/linear_bounds.rs index 1363909..15c628f 100644 --- a/src/linear_bounds.rs +++ b/src/linear_bounds.rs @@ -12,20 +12,20 @@ impl LinearBounds { return LinearBounds {lower, higher}; } - pub fn get_length(self) -> f64 { + pub fn get_length(&self) -> f64 { return (self.higher - self.lower).abs(); } - pub fn get_random_value(self) -> f64 { + pub fn get_random_value(&self) -> f64 { let mut rand: ThreadRng = thread_rng(); return self.lower + self.get_length() * rand.gen::(); } - pub fn get_lower(self) -> f64 { + pub fn get_lower(&self) -> f64 { return self.lower; } - pub fn get_higher(self) -> f64 { + pub fn get_higher(&self) -> f64 { return self.higher; } } \ No newline at end of file diff --git a/src/main.rs b/src/main.rs index cee91a5..b9c474a 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,67 +1,35 @@ -use std::thread::JoinHandle; -use std::time::Instant; use crate::bounds::Bounds; +use crate::function::Function; use crate::linear_bounds::LinearBounds; -use crate::point::Point; - -const SAMPLES: u64 = 1_000_000_000; -const THREAD_CNT: u64 = 8; +use crate::plot::Plot; mod point; mod linear_bounds; mod bounds; +mod function; +mod plot; -fn f(x: f64) -> f64 { - return x.powi(4) - 4_f64 * x.powi(3) + 18_f64 * x.powi(2) - 12_f64 * x - 69_f64 -} - -fn is_inside(point: Point) -> bool { - let y_0: f64 = f(point.get_x()); - return y_0.abs() > point.get_y().abs() && y_0 * point.get_y() >= 0_f64; -} - -fn get_real_integral(bounds: LinearBounds) -> f64 { - return big_f(bounds.get_higher()) - big_f(bounds.get_lower()); -} - -fn big_f(x: f64) -> f64 { - return 0.2 * x.powi(5) - x.powi(4) + 6_f64 * x.powi(3) - 6_f64 * x.powi(2) - 69_f64 * x; -} +const THREAD_CNT: usize = 8; fn main() { - let start: Instant = Instant::now(); - let mut points_inside: u64 = 0; let bounds: Bounds = Bounds::new(LinearBounds::new(0_f64, 20_f64), LinearBounds::new(-100_f64, 150000_f64)); + let func: Function = Function::new(); - let mut threads: Vec> = vec![]; + let sample_limit: usize = 1_000; + let sub_samples: Vec = vec![1, 2, 5]; + let samples_per_iteration: usize = 10; - for _ in 0..THREAD_CNT { - threads.push(std::thread::spawn(move || { - let mut points_inside: u64 = 0; - for _ in 0..SAMPLES { - let point: Point = bounds.get_random_point(); - if is_inside(point) { - points_inside += 1; - } + let mut result: Plot = Plot::new(); + + let mut sample_cnt: usize = 1; + while sample_cnt <= sample_limit { + for sub_sample in sub_samples.clone() { + let samples: usize = sub_sample * sample_cnt; + for _ in 0..samples_per_iteration { + result.add_sample(samples, func.approximate(&bounds, samples, THREAD_CNT)); } - return points_inside; - })) + } + sample_cnt *= 10; } - - for thread in threads { - points_inside += thread.join().unwrap(); - } - - let total_samples: u64 = SAMPLES * THREAD_CNT; - - let approximated_integral: f64 = (points_inside as f64 / total_samples as f64) * bounds.get_area(); - println!("The approximated Integral of the function is: {:.2}", approximated_integral); - - let real_integral: f64 = get_real_integral(bounds.get_x()); - println!("The real Integral of the function is: {:.2}", real_integral); - - let error: f64 = (real_integral - approximated_integral).abs(); - println!("That's an error of {:.2} or {:.5}%", error, (error/real_integral) * 100_f64); - - println!("And the whole thing took {:.5} Seconds for {} samples", start.elapsed().as_micros() as f64 / 1000000_f64, total_samples); -} \ No newline at end of file + result.print(); +} diff --git a/src/plot.rs b/src/plot.rs new file mode 100644 index 0000000..82ab5c2 --- /dev/null +++ b/src/plot.rs @@ -0,0 +1,29 @@ +use std::iter::zip; + +#[derive(Debug)] +pub struct Plot { + x: Vec, + y: Vec>, +} + +impl Plot { + pub fn new() -> Plot { + return Plot { x: vec![], y: vec![] } + } + + pub fn add_sample(&mut self, sample: usize, value: f64) -> () { + if let Some(index) = self.x.iter().position(|&x| x == sample) { + self.y[index].push(value); + } + else { + self.x.push(sample); + self.y.push(vec![value]); + } + } + + pub fn print(self) -> () { + for (x, y) in zip(self.x, self.y) { + println!("{}: {:?}", x, y); + } + } +} \ No newline at end of file diff --git a/src/point.rs b/src/point.rs index 70d0857..210a566 100644 --- a/src/point.rs +++ b/src/point.rs @@ -9,11 +9,11 @@ impl Point { return Point {x, y}; } - pub fn get_x(self) -> f64 { + pub fn get_x(&self) -> f64 { return self.x; } - pub fn get_y(self) -> f64 { + pub fn get_y(&self) -> f64 { return self.y; } } \ No newline at end of file