Add rayon multithreading after intensive juckeling session by Hassan

This commit is contained in:
Jan Bergen 2024-08-02 23:21:16 +02:00
parent 1dbfafc999
commit bd28b6c3e3
4 changed files with 186 additions and 70 deletions

View File

@ -456,6 +456,7 @@ name = "bachelorarbeit"
version = "0.1.0"
dependencies = [
"ahash",
"dashmap",
"eframe",
"egui",
"egui_extras",
@ -465,6 +466,9 @@ dependencies = [
"indicatif",
"nohash-hasher",
"num",
"ray",
"rayon",
"scc",
"wasm-bindgen-futures",
]
@ -849,6 +853,25 @@ dependencies = [
"cfg-if",
]
[[package]]
name = "crossbeam-deque"
version = "0.8.5"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "613f8cc01fe9cf1a3eb3d7f488fd2fa8388403e97039e2f73692932e291a770d"
dependencies = [
"crossbeam-epoch",
"crossbeam-utils",
]
[[package]]
name = "crossbeam-epoch"
version = "0.9.18"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e"
dependencies = [
"crossbeam-utils",
]
[[package]]
name = "crossbeam-utils"
version = "0.8.19"
@ -871,6 +894,20 @@ version = "1.1.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "96a6ac251f4a2aca6b3f91340350eab87ae57c3f127ffeb585e92bd336717991"
[[package]]
name = "dashmap"
version = "6.0.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "804c8821570c3f8b70230c2ba75ffa5c0f9a4189b9a432b6656c536712acae28"
dependencies = [
"cfg-if",
"crossbeam-utils",
"hashbrown",
"lock_api",
"once_cell",
"parking_lot_core",
]
[[package]]
name = "derivative"
version = "2.2.0"
@ -1053,6 +1090,12 @@ dependencies = [
"egui",
]
[[package]]
name = "either"
version = "1.13.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0"
[[package]]
name = "emath"
version = "0.27.2"
@ -2560,6 +2603,32 @@ version = "0.2.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3"
[[package]]
name = "ray"
version = "0.2.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "4706e01695c83e94a0a61af525170acefa9061efc4567dc15ef183bea3b328a8"
[[package]]
name = "rayon"
version = "1.10.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "b418a60154510ca1a002a752ca9714984e21e4241e804d32555251faf8b78ffa"
dependencies = [
"either",
"rayon-core",
]
[[package]]
name = "rayon-core"
version = "1.12.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "1465873a3dfdaa8ae7cb14b4383657caab0b3e8a0aa9ae8e04b044854c8dfce2"
dependencies = [
"crossbeam-deque",
"crossbeam-utils",
]
[[package]]
name = "redox_syscall"
version = "0.3.5"
@ -2679,6 +2748,15 @@ dependencies = [
"winapi-util",
]
[[package]]
name = "scc"
version = "2.1.6"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "05ccfb12511cdb770157ace92d7dda771e498445b78f9886e8cdbc5140a4eced"
dependencies = [
"sdd",
]
[[package]]
name = "scoped-tls"
version = "1.0.1"
@ -2704,6 +2782,12 @@ dependencies = [
"tiny-skia",
]
[[package]]
name = "sdd"
version = "2.1.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "177258b64c0faaa9ffd3c65cd3262c2bc7e2588dbbd9c1641d0346145c1bbda8"
[[package]]
name = "serde"
version = "1.0.199"

View File

@ -15,6 +15,10 @@ egui_plot = "0.27.2"
indicatif = "0.17.8"
nohash-hasher = "0.2.0"
ahash = "0.8.11"
ray = "0.2.0"
rayon = "1.10.0"
scc = "2.1.6"
dashmap = "6.0.1"
[target.'cfg(target_arch = "wasm32")'.dependencies]
wasm-bindgen-futures = "0.4"

View File

@ -1,8 +1,11 @@
use core::f64::consts::PI;
use gauss_quad::GaussLegendre;
use std::collections::HashMap;
use std::io;
use std::rc::Rc;
use std::sync::Arc;
use rayon::prelude::*;
// use scc::DashMap;
use dashmap::DashMap;
use std::{fs::File, io::BufReader};
@ -25,48 +28,53 @@ const INF: f64 = 1e6;
#[derive(Default)]
pub struct Cache {
pub gauss_quad_lut: HashMap<u32, (Rc<[f64]>, Rc<[f64]>), nohash_hasher::BuildNoHashHasher<u32>>,
pub delta_lut: HashMap<u64, f64, nohash_hasher::BuildNoHashHasher<u64>>,
pub omnes_lut: HashMap<
pub gauss_quad_lut: DashMap<u32, (Arc<[f64]>, Arc<[f64]>), nohash_hasher::BuildNoHashHasher<u32>>,
pub delta_lut: DashMap<u64, f64, nohash_hasher::BuildNoHashHasher<u64>>,
pub omnes_lut: DashMap<
u32,
HashMap<u64, Complex, nohash_hasher::BuildNoHashHasher<u64>>,
DashMap<u64, Complex, nohash_hasher::BuildNoHashHasher<u64>>,
nohash_hasher::BuildNoHashHasher<u32>,
>, // First key is n, second is s
}
// pub struct Cache {
// pub gauss_quad_lut: HashMap<u32, (Rc<[f64]>, Rc<[f64]>), RandomState>,
// pub delta_lut: HashMap<u64, f64, RandomState>,
// pub omnes_lut: HashMap<
// pub gauss_quad_lut: DashMap<u32, (Arc<[f64]>, Arc<[f64]>), RandomState>,
// pub delta_lut: DashMap<u64, f64, RandomState>,
// pub omnes_lut: DashMap<
// u32,
// HashMap<u64, Complex, RandomState>,
// DashMap<u64, Complex, RandomState>,
// RandomState,
// >, // First key is n, second is s
// }
// pub struct Cache {
// pub gauss_quad_lut: HashMap<u32, (Rc<[f64]>, Rc<[f64]>)>,
// pub delta_lut: HashMap<u64, f64>,
// pub omnes_lut: HashMap<u32, HashMap<u64, Complex>> // First key is n, second is s
// pub gauss_quad_lut: DashMap<u32, (Arc<[f64]>, Arc<[f64]>)>,
// pub delta_lut: DashMap<u64, f64>,
// pub omnes_lut: DashMap<u32, DashMap<u64, Complex>> // First key is n, second is s
// }
impl Cache {
pub fn from_file(filepath: &str) -> io::Result<Cache> {
// println!("Cache from_file");
let file = File::open(filepath)?;
let mut file = BufReader::new(file);
let me = Cache {
gauss_quad_lut: gq_storage::deserialize(&mut file, |_| true)?
let gql = DashMap::with_hasher(Default::default());
for ent in gq_storage::deserialize(&mut file, |_| true)?
.into_iter()
.map(|(roots, weights)| {
(
u32::try_from(roots.len()).unwrap(),
(Rc::from(roots), Rc::from(weights)),
(Arc::from(roots), Arc::from(weights)),
)
})
.collect(),
delta_lut: HashMap::with_hasher(Default::default()),
omnes_lut: HashMap::with_hasher(Default::default()),
{
let _ = gql.insert(ent.0, ent.1);
}
let me = Cache {
gauss_quad_lut: gql,
delta_lut: DashMap::with_hasher(Default::default()),
omnes_lut: DashMap::with_hasher(Default::default()),
};
Ok(me)
}
@ -74,18 +82,22 @@ impl Cache {
#[allow(unused)]
pub fn from_slice(data: &[u8]) -> io::Result<Cache> {
let mut cursor = std::io::Cursor::new(data);
let me = Cache {
gauss_quad_lut: gq_storage::deserialize(&mut cursor, |_| true)?
let mut gq = DashMap::with_hasher(Default::default());
for ent in gq_storage::deserialize(&mut cursor, |_| true)?
.into_iter()
.map(|(roots, weights)| {
(
u32::try_from(roots.len()).unwrap(),
(Rc::from(roots), Rc::from(weights)),
(Arc::from(roots), Arc::from(weights)),
)
})
.collect(),
delta_lut: HashMap::with_hasher(Default::default()),
omnes_lut: HashMap::with_hasher(Default::default()),
{
gq.insert(ent.0, ent.1);
}
let me = Cache {
gauss_quad_lut: gq,
delta_lut: DashMap::with_hasher(Default::default()),
omnes_lut: DashMap::with_hasher(Default::default()),
};
Ok(me)
}
@ -110,33 +122,32 @@ fn atan_shift(x: f64) -> f64 {
}
}
fn integrate<G: FnMut(f64) -> f64>(
fn integrate<G: Fn(f64) -> f64 + Send + Sync>(
roots: &[f64],
weights: &[f64],
a: f64,
b: f64,
mut f: G,
f: G,
print_prog: bool,
) -> f64 {
let mut sum: f64 = 0.0;
let mut i: usize = 0;
if roots.len() != weights.len() {
panic!("Error: roots and weights are of different length");
}
if print_prog {
println!("Calculating single integration: ");
let prog = indicatif::ProgressBar::new(roots.len() as u64);
while i < roots.len() {
sum += weights[i] * f(roots[i] * (b - a) / 2.0 + (a + b) / 2.0);
i += 1;
println!("Calculating single integration: ");
let prog = indicatif::ProgressBar::new(roots.len() as u64);
let prog_ = prog.clone();
let sum:f64 = (weights.par_iter(), roots.par_iter()).into_par_iter().map(move |(weight,root)| {
let prog = &prog_;
// println!("Orello");
let res = weight * f(root * (b - a) / 2.0 + (a + b) / 2.0);
if print_prog {
prog.inc(1);
}
res
// 0.0
}).sum();
if print_prog {
prog.finish();
} else {
while i < roots.len() {
sum += weights[i] * f(roots[i] * (b - a) / 2.0 + (a + b) / 2.0);
i += 1;
}
}
sum * (b - a) / 2.0
}
@ -148,6 +159,7 @@ fn integrate_complex<G: FnMut(f64) -> Complex>(
b: f64,
mut f: G,
) -> Complex {
// println!("integrate_complex");
let mut sum: Complex = Complex::new(0.0, 0.0);
let mut i: usize = 0;
if roots.len() != weights.len() {
@ -170,12 +182,13 @@ fn k(s: f64) -> f64 {
(s / 4.0 - M_PI.powi(2)).sqrt()
}
pub fn delta_with_lut(cache: &mut Cache, s: f64) -> f64 {
pub fn delta_with_lut(cache: &Cache, s: f64) -> f64 {
// println!("delta_with_lut");
match cache.delta_lut.get(&s.to_bits()) {
Some(val) => *val,
None => {
let val = delta(s);
cache.delta_lut.insert(s.to_bits(), val);
let _ = cache.delta_lut.insert(s.to_bits(), val);
val
}
}
@ -198,51 +211,58 @@ pub fn delta(s: f64) -> f64 {
}
}
pub fn omnes_integrand_tan(cache: &mut Cache, s_tick: f64, s: f64) -> Complex {
pub fn omnes_integrand_tan(cache: &Cache, s_tick: f64, s: f64) -> Complex {
let sub = s_tick.tan();
(delta_with_lut(cache, sub) - delta_with_lut(cache, s))
/ (sub * (sub - s + Complex::new(0.0, EPSILON)) * s_tick.cos().powi(2))
}
pub fn omnes_integrand_notan(cache: &mut Cache, s_tick: f64, s: f64) -> Complex {
pub fn omnes_integrand_notan(cache: &Cache, s_tick: f64, s: f64) -> Complex {
// println!("omnes_integrand_notan");
(delta_with_lut(cache, s_tick) - delta_with_lut(cache, s))
/ (s_tick * (s_tick - s + Complex::new(0.0, EPSILON)))
}
pub fn omnes(cache: &mut Cache, s: f64, n: u32, use_tan: bool) -> Complex {
pub fn omnes(cache: &Cache, s: f64, n: u32, use_tan: bool) -> Complex {
// println!("Omnes");
if let Some(inner_lut) = cache.omnes_lut.get(&n) {
if let Some(res) = inner_lut.get(&s.to_bits()) {
return *res;
}
}
if !cache.omnes_lut.contains_key(&n) {
cache
let _ = cache
.omnes_lut
.insert(n, HashMap::with_hasher(Default::default()));
.insert(n, DashMap::with_hasher(Default::default()));
}
let roots: Rc<[f64]>;
let weights: Rc<[f64]>;
match cache.gauss_quad_lut.get_mut(&n) {
// println!("Retrieving roots and weights omnes");
let roots: Arc<[f64]>;
let weights: Arc<[f64]>;
match cache.gauss_quad_lut.get(&n) {
Some(gq_values) => {
roots = Rc::clone(&gq_values.0);
weights = Rc::clone(&gq_values.1);
// println!("Somobomo");
roots = Arc::clone(&gq_values.0);
weights = Arc::clone(&gq_values.1);
}
None => {
// println!("Nononono");
let gq_values = GaussLegendre::nodes_and_weights(n.try_into().unwrap());
roots = Rc::from(gq_values.0.as_slice());
weights = Rc::from(gq_values.1.as_slice());
roots = Arc::from(gq_values.0.as_slice());
weights = Arc::from(gq_values.1.as_slice());
}
}
// println!("will maybe calc n_omnes");
if !cache.gauss_quad_lut.contains_key(&n) {
let gq_values = GaussLegendre::nodes_and_weights(n.try_into().unwrap());
cache
.gauss_quad_lut
.insert(n, (Rc::from(gq_values.0), Rc::from(gq_values.1)));
println!("n_omnes={} is not included in the lookup file. Consider generating a new lookup file via the appropriate function in gq_storage.rs", n);
let gq_values = GaussLegendre::nodes_and_weights(n.try_into().unwrap());
let _ = cache
.gauss_quad_lut
.insert(n, (Arc::from(gq_values.0), Arc::from(gq_values.1)));
}
// println!("Will calculate integrate complex");
let intgrl = if use_tan {
integrate_complex(&roots, &weights, S0.atan(), PI / 2.0, |s_tick| {
omnes_integrand_tan(cache, s_tick, s)
@ -254,9 +274,9 @@ pub fn omnes(cache: &mut Cache, s: f64, n: u32, use_tan: bool) -> Complex {
};
let omnes_res = (s / PI * (intgrl)).exp()
* (S0 / Complex::new(S0 - s, -EPSILON)).powf(delta_with_lut(cache, s) / PI);
cache
let _ = cache
.omnes_lut
.get_mut(&n)
.get(&n)
.unwrap()
.insert(s.to_bits(), omnes_res);
@ -264,7 +284,7 @@ pub fn omnes(cache: &mut Cache, s: f64, n: u32, use_tan: bool) -> Complex {
}
pub fn phi0_integrand(
cache: &mut Cache,
cache: &Cache,
s_tick: f64,
s: f64,
n_omnes: u32,
@ -277,6 +297,7 @@ pub fn phi0_integrand(
cutoff_s: f64,
cutoff_factor: Complex,
) -> f64 {
// println!("phi0_integrand");
let sub = match (use_tan_phi0, use_xsub) {
(true, true) => s_tick.tan().powi(2) + S0,
(true, false) => s_tick.tan(),
@ -306,7 +327,7 @@ pub fn phi0_integrand(
}
pub fn phi0(
cache: &mut Cache,
cache: &Cache,
s: f64,
n_omnes: u32,
n_phi0: u32,
@ -319,16 +340,20 @@ pub fn phi0(
cutoff_power: f64,
cutoff_s: f64
) -> f64 {
// println!("phi0");
if !cache.gauss_quad_lut.contains_key(&n_phi0) {
let gq_values = GaussLegendre::nodes_and_weights(n_phi0.try_into().unwrap());
cache
let _ = cache
.gauss_quad_lut
.insert(n_phi0, (Rc::from(gq_values.0), Rc::from(gq_values.1)));
.insert(n_phi0, (Arc::from(gq_values.0), Arc::from(gq_values.1)));
println!("n_omnes={} is not included in the lookup file. Consider generating a new lookup file via the appropriate function in gq_storage.rs", n_phi0);
}
let (ref roots, ref weights) = cache.gauss_quad_lut.get(&n_phi0).unwrap();
let (roots, weights) = (Rc::clone(roots), Rc::clone(weights));
let (roots, weights) = {
let ent = cache.gauss_quad_lut.get(&n_phi0).unwrap();
(Arc::clone(&ent.0), Arc::clone(&ent.1))
};
// println!("Will run omnes");
let analyt = if use_reduced_integrand {
omnes(cache, s, n_omnes, use_tan_omnes).norm().ln() / S0.sqrt()
} else {

View File

@ -3,11 +3,13 @@
// - Reduziertes Integral aus Paper zu numerisch berechenbarer Form umformen
use indicatif;
use std::{collections::HashMap, time::Instant};
use std::time::Instant;
use egui_plot::{log_grid_spacer, AxisHints, GridInput, GridMark, Legend, Line, Plot};
use num::{complex::ComplexFloat, Float};
// use scc::HashMap;
mod calc;
mod gq_storage;
@ -129,7 +131,7 @@ impl eframe::App for App {
let omnes_tan_cb =
ui.checkbox(&mut self.use_tan_omnes, "Use tan()-subst. for Omnes ('ot')");
if omnes_tan_cb.clicked() {
self.calc_cache.omnes_lut = HashMap::with_hasher(Default::default()); // TODO split omnes_lut so that this is no longer necessary
self.calc_cache.omnes_lut = dashmap::DashMap::with_hasher(Default::default()); // TODO split omnes_lut so that this is no longer necessary
}
ui.checkbox(&mut self.use_tan_phi0, "Use tan()-subst. for phi0 ('pt')");
ui.checkbox(&mut self.use_xsub_phi0, "Use second subst. for phi0 ('ss')");
@ -595,8 +597,9 @@ fn main() {
// })
// .collect();
let phi0 = calc::phi0(&mut app.calc_cache, 10000.0, 10000, 10000, true, true, true, true, true, false, 0.0, 0.0);
println!("relative difference: {}", phi0 / calc::delta(10000.0) - 1.0);
// let phi0 = calc::phi0(&mut app.calc_cache, 10000.0, 10000, 10000, true, true, true, true, true, false, 0.0, 0.0);
let phi0 = calc::phi0(&mut app.calc_cache, 10000.0, 100000, 100000, true, true, true, true, true, false, 0.0, 0.0);
println!("relative difference: {:e}", phi0 / calc::delta(10000.0) - 1.0);
// eframe::run_native(
// "Omnes Calculator",