From cc61ff3ce4b04577614a121a2ad17e363574f757 Mon Sep 17 00:00:00 2001
From: imrn99 <95699343+imrn99@users.noreply.github.com>
Date: Tue, 15 Oct 2024 17:20:06 +0200
Subject: [PATCH] finish montecarlo impls

---
 integraal/src/structure/implementations.rs | 41 +++++++++++++++++-----
 1 file changed, 32 insertions(+), 9 deletions(-)

diff --git a/integraal/src/structure/implementations.rs b/integraal/src/structure/implementations.rs
index a0ee7f7..93dfd67 100644
--- a/integraal/src/structure/implementations.rs
+++ b/integraal/src/structure/implementations.rs
@@ -5,7 +5,6 @@
 use crate::{
     ComputeMethod, DomainDescriptor, FunctionDescriptor, Integraal, IntegraalError, Scalar,
 };
-use rand::seq::IndexedRandom;
 #[cfg(feature = "montecarlo")]
 use rand::Rng;
 
@@ -161,8 +160,34 @@ fn values_explicit_arm<X: Scalar>(
             ));
         }
         #[cfg(feature = "montecarlo")]
-        ComputeMethod::MonteCarlo { n_sample: _ } => {
-            todo!()
+        ComputeMethod::MonteCarlo { n_sample } => {
+            let (mut min, mut max) = (vals[0], vals[0]);
+            for v in vals {
+                min = min.min(*v);
+                max = max.max(*v);
+            }
+            let height: X = max - min;
+            let widths = args.windows(2).enumerate().map(|(i, slice)| {
+                let [a, b] = slice else { unreachable!() };
+                (vals[i].min(X::zero())..vals[i].max(X::zero()), (*b - *a))
+            });
+            let mut rng = rand::thread_rng();
+            let random_numbers: Vec<X> = (&mut rng)
+                .sample_iter(
+                    rand::distr::Uniform::new(min.to_f64().unwrap(), max.to_f64().unwrap())
+                        .unwrap(),
+                )
+                .take(*n_sample * (args.len() - 1))
+                .map(|s| X::from(s).unwrap())
+                .collect();
+            let total_in: X = random_numbers
+                .chunks_exact(*n_sample)
+                .zip(widths)
+                .map(|(samples, (range, width))| {
+                    X::from(samples.iter().filter(|s| range.contains(s)).count()).unwrap() * width
+                })
+                .sum();
+            height * total_in / X::from(*n_sample as f64).unwrap()
         }
     };
 
@@ -370,12 +395,10 @@ fn closure_explicit_arm<X: Scalar>(
         }
         #[cfg(feature = "montecarlo")]
         ComputeMethod::MonteCarlo { n_sample } => {
-            // FIXME: nuke this temp allocation
-            let vals: Vec<_> = args[..args.len() - 2].iter().map(|x| closure(*x)).collect();
-            let (mut min, mut max) = (vals[0], vals[0]);
-            for v in &vals {
-                min = min.min(*v);
-                max = max.max(*v);
+            let (mut min, mut max) = (closure(args[0]), closure(args[0]));
+            for a in args {
+                min = min.min(closure(*a));
+                max = max.max(closure(*a));
             }
             let height: X = max - min;
             let widths = args.windows(2).map(|slice| {