Skip to content

Commit

Permalink
feat: implement Monte-Carlo compute method (#36)
Browse files Browse the repository at this point in the history 
* first attempt at montecarlo

* implement montecarlo for closure/uniform inputs

* implement motnecarlo for closure/explicit inputs

* finish montecarlo impls

* fix missing division

* add tests for montecarlo compute
  • Loading branch information
@imrn99
imrn99 authored Oct 15, 2024
1 parent dfc0843 commit e0a8d26
Show file tree
Hide file tree
Showing 7 changed files with 170 additions and 9 deletions.
121 changes: 112 additions & 9 deletions integraal/src/structure/implementations.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@
use crate::{
ComputeMethod, DomainDescriptor, FunctionDescriptor, Integraal, IntegraalError, Scalar,
};
#[cfg(feature = "montecarlo")]
use rand::Rng;

// ------ CONTENT

Expand Down Expand Up @@ -158,15 +160,41 @@ 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 * (args.len() - 1)).unwrap()
}
};

Ok(res)
}

#[allow(clippy::cast_possible_truncation)]
#[allow(clippy::cast_possible_truncation, clippy::too_many_lines)]
fn values_uniform_arm<X: Scalar>(
vals: &[X],
domain: &DomainDescriptor<X>,
Expand Down Expand Up @@ -278,8 +306,29 @@ fn values_uniform_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 intervals = vals.iter().map(|v| v.min(X::zero())..v.max(X::zero()));
let volume: X = (max - min) * (X::from(*n_step).unwrap() * *step);
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 * *n_step)
.map(|s| X::from(s).unwrap())
.collect();
let total_in: usize = random_numbers
.chunks_exact(*n_sample)
.zip(intervals)
.map(|(samples, range)| samples.iter().filter(|s| range.contains(s)).count())
.sum();
volume * X::from(total_in as f64 / (*n_sample * *n_step) as f64).unwrap()
}
};

Expand Down Expand Up @@ -345,8 +394,37 @@ fn closure_explicit_arm<X: Scalar>(
));
}
#[cfg(feature = "montecarlo")]
ComputeMethod::MonteCarlo { n_sample: _ } => {
todo!()
ComputeMethod::MonteCarlo { n_sample } => {
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| {
let [a, b] = slice else { unreachable!() };
(
closure(*a).min(X::zero())..closure(*a).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 * (args.len() - 1)).unwrap()
}
};
Ok(res)
Expand Down Expand Up @@ -476,8 +554,33 @@ fn closure_uniform_arm<X: Scalar>(
}
}
#[cfg(feature = "montecarlo")]
ComputeMethod::MonteCarlo { n_sample: _ } => {
todo!()
ComputeMethod::MonteCarlo { n_sample } => {
// FIXME: nuke this temp allocation
let vals: Vec<_> = (0..*n_step)
.map(|i| closure(*start + X::from(i).unwrap() * *step))
.collect();
let (mut min, mut max) = (vals[0], vals[0]);
for v in &vals {
min = min.min(*v);
max = max.max(*v);
}
let intervals = vals.iter().map(|v| v.min(X::zero())..v.max(X::zero()));
let volume: X = (max - min) * (X::from(*n_step).unwrap() * *step);
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 * *n_step)
.map(|s| X::from(s).unwrap())
.collect();
let total_in: usize = random_numbers
.chunks_exact(*n_sample)
.zip(intervals)
.map(|(samples, range)| samples.iter().filter(|s| range.contains(s)).count())
.sum();
volume * X::from(total_in as f64 / (*n_sample * *n_step) as f64).unwrap()
}
};

Expand Down
2 changes: 2 additions & 0 deletions integraal/src/tests/function_a.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ mod double {
const SIMPSON_TOLERANCE: f64 = 1e-5;
const BOOLE_TOLERANCE: f64 = 1e-5;
const ROMBERG_TOLERANCE: f64 = 1e-5;
const MONTECARLO_TOLERANCE: f64 = 1e-5;

all_tests!(
f64,
Expand Down Expand Up @@ -50,6 +51,7 @@ mod simple {
const SIMPSON_TOLERANCE: f32 = 1e-5;
const BOOLE_TOLERANCE: f32 = 1e-5;
const ROMBERG_TOLERANCE: f32 = 1e-5;
const MONTECARLO_TOLERANCE: f32 = 1e-5;

all_tests!(
f32,
Expand Down
2 changes: 2 additions & 0 deletions integraal/src/tests/function_b.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ mod double {
const SIMPSON_TOLERANCE: f64 = 1e-5;
const BOOLE_TOLERANCE: f64 = 1e-5;
const ROMBERG_TOLERANCE: f64 = 1e-5;
const MONTECARLO_TOLERANCE: f64 = 1e-5;

all_tests!(
f64,
Expand Down Expand Up @@ -64,6 +65,7 @@ mod simple {
const SIMPSON_TOLERANCE: f32 = 1e-5;
const BOOLE_TOLERANCE: f32 = 1e-5;
const ROMBERG_TOLERANCE: f32 = 1e-5;
const MONTECARLO_TOLERANCE: f32 = 1e-5;

all_tests!(
f32,
Expand Down
2 changes: 2 additions & 0 deletions integraal/src/tests/function_c.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ mod double {
const SIMPSON_TOLERANCE: f64 = 1e-5;
const BOOLE_TOLERANCE: f64 = 1e-5;
const ROMBERG_TOLERANCE: f64 = 1e-5;
const MONTECARLO_TOLERANCE: f64 = 1e-5;

all_tests!(
f64,
Expand Down Expand Up @@ -62,6 +63,7 @@ mod simple {
const SIMPSON_TOLERANCE: f32 = 1e-5;
const BOOLE_TOLERANCE: f32 = 1e-5;
const ROMBERG_TOLERANCE: f32 = 1e-5;
const MONTECARLO_TOLERANCE: f32 = 1e-5;

all_tests!(
f32,
Expand Down
2 changes: 2 additions & 0 deletions integraal/src/tests/function_d.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ mod double {
const SIMPSON_TOLERANCE: f64 = 1e-5;
const BOOLE_TOLERANCE: f64 = 1e-5;
const ROMBERG_TOLERANCE: f64 = 1e-5;
const MONTECARLO_TOLERANCE: f64 = 1e-5;

all_tests!(
f64,
Expand Down Expand Up @@ -51,6 +52,7 @@ mod simple {
const SIMPSON_TOLERANCE: f32 = 1e-5;
const BOOLE_TOLERANCE: f32 = 1e-5;
const ROMBERG_TOLERANCE: f32 = 1e-5;
const MONTECARLO_TOLERANCE: f32 = 1e-5;

all_tests!(
f32,
Expand Down
2 changes: 2 additions & 0 deletions integraal/src/tests/function_e.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ mod double {
const SIMPSON_TOLERANCE: f64 = 1e-5;
const BOOLE_TOLERANCE: f64 = 1e-5;
const ROMBERG_TOLERANCE: f64 = 1e-5;
const MONTECARLO_TOLERANCE: f64 = 1e-5;

all_tests!(
f64,
Expand Down Expand Up @@ -51,6 +52,7 @@ mod simple {
const SIMPSON_TOLERANCE: f32 = 1e-5;
const BOOLE_TOLERANCE: f32 = 1e-5;
const ROMBERG_TOLERANCE: f32 = 1e-5;
const MONTECARLO_TOLERANCE: f32 = 1e-5;

all_tests!(
f32,
Expand Down
48 changes: 48 additions & 0 deletions integraal/src/tests/mod.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -360,6 +360,54 @@ macro_rules! all_tests {
ROMBERG_TOLERANCE
);
}

#[cfg(feature = "montecarlo")]
#[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
mod montecarlo {
use super::*;

generate_test!(
$ft,
ClosureExplicit,
$dm,
$fnd_cls,
$dmd_xpl,
ComputeMethod::MonteCarlo { n_sample: 100 },
RES,
MONTECARLO_TOLERANCE
);

generate_test!(
$ft,
ClosureUniform,
$fnd_cls,
$dmd_uni,
ComputeMethod::MonteCarlo { n_sample: 100 },
RES,
MONTECARLO_TOLERANCE
);

generate_test!(
$ft,
ValuesExplicit,
$dm,
$fnd_val,
$dmd_xpl,
ComputeMethod::MonteCarlo { n_sample: 100 },
RES,
MONTECARLO_TOLERANCE
);

generate_test!(
$ft,
ValuesUniform,
$fnd_val,
$dmd_uni,
ComputeMethod::MonteCarlo { n_sample: 100 },
RES,
MONTECARLO_TOLERANCE
);
}
};
}

Expand Down

0 comments on commit e0a8d26

Please sign in to comment.