jan/bachelor
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Push current version

Browse Source
This commit is contained in:
Jan Bergen 2024-07-16 15:48:35 +02:00
parent c1dd180c61
commit f520c8b3cb
3 changed files with 519 additions and 271 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

313
bachelorarbeit/src/calc.rs
View File

@ -1,21 +1,16 @@
use core::f64::consts::PI; use core::f64::consts::PI;
use std::time::Instant;
use gauss_quad::GaussLegendre; use gauss_quad::GaussLegendre;
use num::Float;
use std::collections::HashMap; use std::collections::HashMap;
use std::io; use std::io;
use std::rc::Rc; use std::rc::Rc;
use std::{ use std::{fs::File, io::BufReader};
fs::File,
io::BufReader,
};
use crate::gq_storage; // This should be moved to main.rs after new archtitecture (see main.rs TODO) use crate::gq_storage;
pub type Complex = num::complex::Complex<f64>; pub type Complex = num::complex::Complex<f64>;
const M_ROH: f64 = 0.7736; const M_ROH: f64 = 0.7736;
const M_PI: f64 = 0.13957; const M_PI: f64 = 0.13957;
pub const S0: f64 = 4.0 * M_PI * M_PI; pub const S0: f64 = 4.0 * M_PI * M_PI;
const M_K: f64 = 0.496; const M_K: f64 = 0.496;
@ -26,35 +21,55 @@ const LAMBDA2: f64 = -1.96;
const LAMBDA_HIGH: f64 = 10.0; const LAMBDA_HIGH: f64 = 10.0;
pub const EPSILON: f64 = 1e-15; pub const EPSILON: f64 = 1e-15;
const S_MID_CUTOFF: f64 = 1.42 * 1.42; const S_MID_CUTOFF: f64 = 1.42 * 1.42;
const INF: f64 = 1e5; const INF: f64 = 1e6;
#[derive(Default)] #[derive(Default)]
pub struct Cache { pub struct Cache {
pub gauss_quad_lut: HashMap<u32, (Rc<[f64]>, Rc<[f64]>), nohash_hasher::BuildNoHashHasher<u32>>, 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 delta_lut: HashMap<u64, f64, nohash_hasher::BuildNoHashHasher<u64>>,
pub omnes_lut: HashMap<u32, HashMap<u64, Complex, nohash_hasher::BuildNoHashHasher<u64>>, nohash_hasher::BuildNoHashHasher<u32>> // First key is n, second is s pub omnes_lut: HashMap<
u32,
HashMap<u64, Complex, nohash_hasher::BuildNoHashHasher<u64>>,
nohash_hasher::BuildNoHashHasher<u32>,
>, // First key is n, second is s
} }
impl Cache { impl Cache {
pub fn from_file(filepath: &str) -> io::Result<Cache> { pub fn from_file(filepath: &str) -> io::Result<Cache> {
let file = File::open(filepath)?; let file = File::open(filepath)?;
let mut file = BufReader::new(file); let mut file = BufReader::new(file);
let me = Cache{ let me = Cache {
gauss_quad_lut: gq_storage::deserialize(&mut file, |_| true)?.into_iter() gauss_quad_lut: gq_storage::deserialize(&mut file, |_| true)?
.map(|(roots, weights)| (u32::try_from(roots.len()).unwrap(), (Rc::from(roots), Rc::from(weights)))).collect(), .into_iter()
.map(|(roots, weights)| {
(
u32::try_from(roots.len()).unwrap(),
(Rc::from(roots), Rc::from(weights)),
)
})
.collect(),
delta_lut: HashMap::with_hasher(Default::default()), delta_lut: HashMap::with_hasher(Default::default()),
omnes_lut: HashMap::with_hasher(Default::default()) omnes_lut: HashMap::with_hasher(Default::default()),
}; };
Ok(me) Ok(me)
} }
#[allow(unused)]
pub fn from_slice(data: &[u8]) -> io::Result<Cache> { pub fn from_slice(data: &[u8]) -> io::Result<Cache> {
let mut cursor = std::io::Cursor::new(data); let mut cursor = std::io::Cursor::new(data);
let me = Cache{ let me = Cache {
gauss_quad_lut: gq_storage::deserialize(&mut cursor, |_| true)?.into_iter() gauss_quad_lut: gq_storage::deserialize(&mut cursor, |_| true)?
.map(|(roots, weights)| (u32::try_from(roots.len()).unwrap(), (Rc::from(roots), Rc::from(weights)))).collect(), .into_iter()
.map(|(roots, weights)| {
(
u32::try_from(roots.len()).unwrap(),
(Rc::from(roots), Rc::from(weights)),
)
})
.collect(),
delta_lut: HashMap::with_hasher(Default::default()), delta_lut: HashMap::with_hasher(Default::default()),
omnes_lut: HashMap::with_hasher(Default::default()) omnes_lut: HashMap::with_hasher(Default::default()),
}; };
Ok(me) Ok(me)
} }
@ -63,7 +78,7 @@ impl Cache {
// derivated values // derivated values
fn lambda_0() -> f64 { fn lambda_0() -> f64 {
delta((2.0*M_K).powi(2)) delta((2.0 * M_K).powi(2))
} }
fn delta_mid_cutoff() -> f64 { fn delta_mid_cutoff() -> f64 {
@ -79,31 +94,54 @@ fn atan_shift(x: f64) -> f64 {
} }
} }
fn integrate<G:FnMut(f64) -> f64>(roots: &[f64], weights: &[f64], a:f64, b:f64, mut f: G) -> f64 { fn integrate<G: FnMut(f64) -> f64>(
roots: &[f64],
weights: &[f64],
a: f64,
b: f64,
mut f: G,
print_prog: bool,
) -> f64 {
let mut sum: f64 = 0.0; let mut sum: f64 = 0.0;
let mut i: usize = 0; let mut i: usize = 0;
if roots.len() != weights.len() { if roots.len() != weights.len() {
panic!("Error: roots and weights are of different length"); panic!("Error: roots and weights are of different length");
} }
if print_prog {
while i < roots.len() { println!("Calculating single integration: ");
sum += weights[i] * f(roots[i]*(b-a)/2.0 + (a+b)/2.0); let prog = indicatif::ProgressBar::new(roots.len() as u64);
i += 1; while i < roots.len() {
sum += weights[i] * f(roots[i] * (b - a) / 2.0 + (a + b) / 2.0);
i += 1;
prog.inc(1);
}
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 sum * (b - a) / 2.0
} }
fn integrate_complex<G:FnMut(f64) -> Complex>(roots: &[f64], weights: &[f64], a:f64, b:f64, mut f: G) -> Complex { fn integrate_complex<G: FnMut(f64) -> Complex>(
let mut sum: Complex = Complex::new(0.0,0.0); roots: &[f64],
weights: &[f64],
a: f64,
b: f64,
mut f: G,
) -> Complex {
let mut sum: Complex = Complex::new(0.0, 0.0);
let mut i: usize = 0; let mut i: usize = 0;
if roots.len() != weights.len() { if roots.len() != weights.len() {
panic!("Error: roots and weights are of different length"); panic!("Error: roots and weights are of different length");
} }
while i < roots.len() { while i < roots.len() {
sum += weights[i] * f(roots[i]*(b-a)/2.0 + (a+b)/2.0); sum += weights[i] * f(roots[i] * (b - a) / 2.0 + (a + b) / 2.0);
i += 1; i += 1;
} }
sum * (b-a)/2.0 sum * (b - a) / 2.0
} }
fn omega(s: f64) -> f64 { fn omega(s: f64) -> f64 {
@ -113,14 +151,12 @@ fn omega(s: f64) -> f64 {
} }
fn k(s: f64) -> f64 { fn k(s: f64) -> f64 {
(s/4.0 - M_PI.powi(2)).sqrt() (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: &mut Cache, s: f64) -> f64 {
match cache.delta_lut.get(&s.to_bits()) { match cache.delta_lut.get(&s.to_bits()) {
Some(val) => { Some(val) => *val,
*val
}
None => { None => {
let val = delta(s); let val = delta(s);
cache.delta_lut.insert(s.to_bits(), val); cache.delta_lut.insert(s.to_bits(), val);
@ -132,154 +168,179 @@ pub fn delta_with_lut(cache:&mut Cache, s: f64) -> f64 {
pub fn delta(s: f64) -> f64 { pub fn delta(s: f64) -> f64 {
if s <= 2.0 * M_K { if s <= 2.0 * M_K {
atan_shift( atan_shift(
(s.sqrt() / (2.0*k(s).powi(3)) (s.sqrt() / (2.0 * k(s).powi(3))
* (M_ROH.powi(2) - s) * (M_ROH.powi(2) - s)
* ( (2.0*M_PI.powi(3)) / (M_ROH.powi(2) * s.sqrt()) + B0 + B1 * omega(s)) * ((2.0 * M_PI.powi(3)) / (M_ROH.powi(2) * s.sqrt()) + B0 + B1 * omega(s)))
).powi(-1) .powi(-1),
) )
} else if s <= S_MID_CUTOFF { } else if s <= S_MID_CUTOFF {
let par = 0.5 * s.sqrt() / M_K - 1.0; let par = 0.5 * s.sqrt() / M_K - 1.0;
lambda_0() + LAMBDA1 * par + LAMBDA2 * par.powi(2) lambda_0() + LAMBDA1 * par + LAMBDA2 * par.powi(2)
} else { } else {
PI + PI + (delta_mid_cutoff() - PI) * (LAMBDA_HIGH * LAMBDA_HIGH + S_MID_CUTOFF)
(delta_mid_cutoff() - PI) / (LAMBDA_HIGH * LAMBDA_HIGH + s)
* (LAMBDA_HIGH * LAMBDA_HIGH + S_MID_CUTOFF)
/ (LAMBDA_HIGH * LAMBDA_HIGH + s)
} }
} }
pub fn omnes_integrand_tan(cache:&mut Cache, s_tick:f64, s:f64) -> Complex { pub fn omnes_integrand_tan(cache: &mut Cache, s_tick: f64, s: f64) -> Complex {
let sub = s_tick.tan(); let sub = s_tick.tan();
(delta_with_lut(cache, sub) - delta_with_lut(cache, s)) (delta_with_lut(cache, sub) - delta_with_lut(cache, s))
/ ( / (sub * (sub - s + Complex::new(0.0, EPSILON)) * s_tick.cos().powi(2))
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: &mut Cache, s_tick: f64, s: f64) -> Complex {
(delta_with_lut(cache, s_tick) - delta_with_lut(cache, s)) (delta_with_lut(cache, s_tick) - delta_with_lut(cache, s))
/ (s_tick * (s_tick - s + Complex::new(0.0, EPSILON))) / (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: &mut Cache, s: f64, n: u32, use_tan: bool) -> Complex {
let roots: Rc<[f64]>; let roots: Rc<[f64]>;
let weights: Rc<[f64]>; let weights: Rc<[f64]>;
if let Some(inner_lut) = cache.omnes_lut.get(&n) { if let Some(inner_lut) = cache.omnes_lut.get(&n) {
if let Some(res) = inner_lut.get(&s.to_bits()) { if let Some(res) = inner_lut.get(&s.to_bits()) {
return *res return *res;
} }
} }
if !cache.omnes_lut.contains_key(&n) { if !cache.omnes_lut.contains_key(&n) {
cache.omnes_lut.insert(n, HashMap::with_hasher(Default::default())); cache
} .omnes_lut
.insert(n, HashMap::with_hasher(Default::default()));
}
match cache.gauss_quad_lut.get_mut(&n) { match cache.gauss_quad_lut.get_mut(&n) {
Some(gq_values) => { Some(gq_values) => {
roots = Rc::clone(&gq_values.0); roots = Rc::clone(&gq_values.0);
weights = Rc::clone(&gq_values.1); weights = Rc::clone(&gq_values.1);
} }
None => { None => {
let gq_values = GaussLegendre::nodes_and_weights(n.try_into().unwrap()); let gq_values = GaussLegendre::nodes_and_weights(n.try_into().unwrap());
roots = Rc::from(gq_values.0.as_slice()); roots = Rc::from(gq_values.0.as_slice());
weights = Rc::from(gq_values.1.as_slice()); weights = Rc::from(gq_values.1.as_slice());
} }
} }
if !cache.gauss_quad_lut.contains_key(&n) { if !cache.gauss_quad_lut.contains_key(&n) {
let gq_values = GaussLegendre::nodes_and_weights(n.try_into().unwrap()); 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))); cache
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); .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 intgrl = if use_tan { let intgrl = if use_tan {
integrate_complex(&roots, &weights, S0.atan(), PI / 2.0, |s_tick| omnes_integrand_tan(cache, s_tick, s)) integrate_complex(&roots, &weights, S0.atan(), PI / 2.0, |s_tick| {
omnes_integrand_tan(cache, s_tick, s)
})
} else { } else {
integrate_complex(&roots, &weights, S0, INF, |s_tick| omnes_integrand_notan(cache, s_tick, s)) integrate_complex(&roots, &weights, S0, INF, |s_tick| {
omnes_integrand_notan(cache, s_tick, s)
})
}; };
let omnes_res = (s/PI*(intgrl)).exp() * (S0 / Complex::new(S0 - s, -EPSILON)).powf(delta_with_lut(cache, s) / PI); let omnes_res = (s / PI * (intgrl)).exp()
cache.omnes_lut.get_mut(&n).unwrap().insert(s.to_bits(), omnes_res); * (S0 / Complex::new(S0 - s, -EPSILON)).powf(delta_with_lut(cache, s) / PI);
cache
.omnes_lut
.get_mut(&n)
.unwrap()
.insert(s.to_bits(), omnes_res);
omnes_res omnes_res
} }
pub fn phi0_integrand_tan_xsub(cache:&mut Cache, s_tick: f64, s: f64, n_omnes:u32, use_tan_omnes: bool) -> f64 { pub fn phi0_integrand(
let x_sub = s_tick.tan().powi(2) + S0; cache: &mut Cache,
s_tick: f64,
s: f64,
n_omnes: u32,
use_tan_omnes: bool,
use_tan_phi0: bool,
use_xsub: bool,
_use_reduced_integrand: bool,
) -> f64 {
let sub = match (use_tan_phi0, use_xsub) {
(true, true) => s_tick.tan().powi(2) + S0,
(true, false) => s_tick.tan(),
(false, true) => s_tick.powi(2) + S0,
(false, false) => s_tick,
};
let omnes_s = omnes(cache, s, n_omnes, use_tan_omnes); let omnes_s = omnes(cache, s, n_omnes, use_tan_omnes);
let omnes_x = omnes(cache, x_sub, n_omnes, use_tan_omnes); let omnes_sub = omnes(cache, sub, n_omnes, use_tan_omnes);
(omnes_x / omnes_s).norm_sqr().ln() let mut res = (omnes_sub / omnes_s).norm().ln() / (sub * (sub - s));
/ ( (x_sub)* (x_sub-s) ) res *= match (use_tan_phi0, use_xsub) {
/ s_tick.cos().powi(2) // Don't use x_sub here, since this term is from tan substitution itself (true, true) => 2.0 / s_tick.cos().powi(2),
} (true, false) => 1.0 / ((sub - S0).sqrt() * s_tick.cos().powi(2)),
(false, true) => 2.0,
pub fn phi0_integrand_notan_xsub(cache:&mut Cache, s_tick: f64, s: f64, n_omnes:u32, use_tan_omnes: bool) -> f64 { (false, false) => 1.0 / (sub - S0).sqrt(),
let x_sub = s_tick.powi(2) + S0; };
let omnes_s = omnes(cache, s, n_omnes, use_tan_omnes);
let omnes_x = omnes(cache, x_sub, n_omnes, use_tan_omnes);
(omnes_x / omnes_s).norm_sqr().ln()
/ ( (x_sub)* (x_sub-s) )
}
pub fn phi0_integrand_tan_noxsub(cache:&mut Cache, s_tick: f64, s: f64, n_omnes:u32, use_tan_omnes: bool) -> f64 {
let tansub = s_tick.tan();
println!("lol {:?}", omnes(cache, tansub, n_omnes, use_tan_omnes));
let omnes_s = omnes(cache, s, n_omnes, use_tan_omnes);
let omnes_s_tick = omnes(cache, tansub, n_omnes, use_tan_omnes);
let res = (omnes_s_tick / omnes_s).norm_sqr().ln()
/ ( (tansub - S0).sqrt() * tansub * (tansub - s) )
/ s_tick.cos().powi(2);
println!("{:?} {:?} {:?} {:?} {:?}", s_tick, tansub, omnes_s, omnes_s_tick, res);
println!("LOL {:?}", omnes(cache, tansub, n_omnes, use_tan_omnes));
res res
} }
pub fn phi0_integrand_notan_noxsub(cache:&mut Cache, s: f64, s_tick: f64, n_omnes:u32, use_tan_omnes: bool) -> f64 { pub fn phi0(
cache: &mut Cache,
let omnes_s = omnes(cache, s, n_omnes, use_tan_omnes); s: f64,
let omnes_s_tick = omnes(cache, s_tick, n_omnes, use_tan_omnes); n_omnes: u32,
let res = (omnes_s_tick / omnes_s).norm_sqr().ln() n_phi0: u32,
/ ( (s_tick - S0).sqrt() * s_tick * (s_tick - s) ); use_tan_omnes: bool,
println!("{:?} {:?} {:?} {:?} {:?}", s, s_tick, omnes_s, omnes_s_tick, res); use_tan_phi0: bool,
use_xsub: bool,
res use_reduced_integrand: bool,
print_prog: bool,
} ) -> f64 {
pub fn phi0(cache:&mut Cache, s: f64, n_omnes: u32, n_phi0: u32, use_tan_omnes: bool, use_tan_phi0: bool, use_xsub: bool) -> f64 {
//let t0 = Instant::now();
if !cache.gauss_quad_lut.contains_key(&n_phi0) { if !cache.gauss_quad_lut.contains_key(&n_phi0) {
let gq_values = GaussLegendre::nodes_and_weights(n_phi0.try_into().unwrap()); let gq_values = GaussLegendre::nodes_and_weights(n_phi0.try_into().unwrap());
cache.gauss_quad_lut.insert(n_phi0, (Rc::from(gq_values.0), Rc::from(gq_values.1))); cache
.gauss_quad_lut
.insert(n_phi0, (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_phi0); 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 (ref roots, ref weights) = cache.gauss_quad_lut.get(&n_phi0).unwrap();
let (roots, weights) = (Rc::clone(roots), Rc::clone(weights)); let (roots, weights) = (Rc::clone(roots), Rc::clone(weights));
let analyt = omnes(cache, s, n_omnes, use_tan_omnes).norm_sqr().ln() / S0.sqrt() / 2.0; let analyt = if use_reduced_integrand {
let intgrl = if use_tan_phi0 { omnes(cache, s, n_omnes, use_tan_omnes).norm().ln() / S0.sqrt()
if use_xsub {
s / PI * integrate(&*roots, &*weights, 0.0, PI / 2.0, |x| phi0_integrand_tan_xsub(cache, x, s, n_omnes, use_tan_omnes))
} else {
s / PI * integrate(&*roots, &*weights, S0.atan(), PI / 2.0, |x| phi0_integrand_tan_noxsub(cache, x, s, n_omnes, use_tan_omnes))
}
} else { } else {
if use_xsub { (omnes(cache, s, n_omnes, use_tan_omnes) / omnes(cache, S0, n_omnes, use_tan_omnes))
s / PI * integrate(&*roots, &*weights, 0.0, INF, |x| phi0_integrand_notan_xsub(cache, x, s, n_omnes, use_tan_omnes)) .norm()
} else { .ln()
s / PI * integrate(&*roots, &*weights, S0, INF, |x| phi0_integrand_notan_noxsub(cache, x, s, n_omnes, use_tan_omnes)) / (s * (s - S0).powf(1.5))
}
}; };
//println!("{:?}", t0.elapsed()); let intgrl = integrate(
-(s-S0).sqrt() * (intgrl - analyt) &*roots,
&*weights,
if use_xsub {
0.0
} else if use_tan_phi0 {
S0.atan()
} else {
S0
},
if use_tan_phi0 { PI / 2.0 } else { INF },
|x| {
phi0_integrand(
cache,
x,
s,
n_omnes,
use_tan_omnes,
use_tan_phi0,
use_xsub,
use_reduced_integrand,
) * s
/ PI
},
print_prog,
);
-(s - S0).sqrt() * (intgrl - analyt)
} }
#[allow(unused)] #[allow(unused)]
pub fn bit_pattern_randomness_tester(a:f64, b:f64, n_points: u64) -> Vec<f64>{ pub fn bit_pattern_randomness_tester(a: f64, b: f64, n_points: u64) -> Vec<f64> {
let mut modulos = vec![0.0; n_points as usize]; let mut modulos = vec![0.0; n_points as usize];
let x_values: Vec<u64> = (0..modulos.len() as u32).map(|x| (f64::from(x) * ((b - a) / f64::from(modulos.len() as u32)) + a).to_bits()).collect(); let x_values: Vec<u64> = (0..modulos.len() as u32)
.map(|x| (f64::from(x) * ((b - a) / f64::from(modulos.len() as u32)) + a).to_bits())
.collect();
for i in 0..x_values.len() { for i in 0..x_values.len() {
let val = (x_values[i] % n_points) as usize; let val = (x_values[i] % n_points) as usize;
modulos[val] += 1.0; modulos[val] += 1.0;

2
bachelorarbeit/src/gq_storage.rs
View File

@ -28,7 +28,7 @@ pub fn serialize(filepath: &str, params: Vec<(Vec<f64>, Vec<f64>)>) -> io::Resul
Ok(()) Ok(())
} }
pub fn deserialize<R:std::io::Read + std::io::Seek>(reader: &mut R, filter: impl Fn(u64) -> bool) -> io::Result<Vec<(Vec<f64>, Vec<f64>)>> { pub fn deserialize<R:Read + Seek>(reader: &mut R, filter: impl Fn(u64) -> bool) -> io::Result<Vec<(Vec<f64>, Vec<f64>)>> {
let mut buf = [0u8; 8]; let mut buf = [0u8; 8];
let mut values = Vec::new(); let mut values = Vec::new();

475
bachelorarbeit/src/main.rs
View File

@ -1,38 +1,35 @@
// TODO // TODO
// - Lookup Table implementieren -> Weights, roots schon früher rausziehen?
// -> Warum macht delta_lut alles langsamer?
// - Omnes auch komplex integrieren (geht nicht mit Integralmethode von gauss-quad)
// - Textfelder und "Berechnen"-Button in UI, um z.B. Wertebereich zu setzen
// - Neue Architketur: struct calculator, das lookups speichert -> In die nötigen Funktionen übergeben. Dadurch kann man omnes auch als single calc machen und trotzdem die Punkte nur einmal berechnen
// - AutoImport ausschalten // - AutoImport ausschalten
// - Reduziertes Integral aus Paper zu numerisch berechenbarer Form umformen
use indicatif;
use std::collections::HashMap;
use egui_plot::{log_grid_spacer, AxisHints, GridInput, GridMark, Legend, Line, Plot}; use egui_plot::{log_grid_spacer, AxisHints, GridInput, GridMark, Legend, Line, Plot};
use num::complex::ComplexFloat; use num::complex::ComplexFloat;
use std::time::Instant;
mod gq_storage;
mod calc; mod calc;
mod gq_storage;
struct App { struct App {
plot_data: Vec<PlotCurve>, plot_data: Vec<PlotCurve>,
plots_available: [String; 6], plots_available: [String; 4],
calc_cache: calc::Cache, calc_cache: calc::Cache,
scaling_type: ScalingType, scaling_type: ScalingType,
a: f64, a: f64,
b: f64, b: f64,
n_points: u32, n_points: u32,
n_omnes: u32, n_omnes: u32,
n_phi0: u32, n_phi0: u32,
use_tan_omnes: bool, use_tan_omnes: bool,
use_tan_phi0: bool, use_tan_phi0: bool,
use_xsub_phi0: bool use_xsub_phi0: bool,
use_reduced_integrand: bool,
} }
#[derive(Clone)] #[derive(Clone)]
struct PlotCurve { struct PlotCurve {
points: Vec<[f64;2]>, points: Vec<[f64; 2]>,
name: String, name: String,
} }
@ -40,26 +37,32 @@ impl Default for App {
fn default() -> Self { fn default() -> Self {
App { App {
plot_data: Vec::default(), plot_data: Vec::default(),
plots_available: ["Delta".to_string(), "Omnes Integrand".to_string(), "Omnes".to_string(), "Phi0 Integrand".to_string(), "Phi0".to_string(), "Phi0 / delta".to_string()], plots_available: [
"Delta".to_string(),
"Omnes".to_string(),
"Phi0".to_string(),
"Phi0 / delta".to_string(),
],
calc_cache: calc::Cache::default(), calc_cache: calc::Cache::default(),
scaling_type: ScalingType::default(), scaling_type: ScalingType::default(),
a: calc::S0, a: calc::S0,
b: 9.0, b: 9.0,
n_points: 5000, n_points: 500,
n_omnes:1000, n_omnes: 500,
n_phi0:1000, n_phi0: 500,
use_tan_omnes: true, use_tan_omnes: true,
use_tan_phi0: true, use_tan_phi0: true,
use_xsub_phi0: true use_xsub_phi0: true,
use_reduced_integrand: true,
} }
} }
} }
#[derive(Debug, Default, PartialEq, Clone, Copy)] #[derive(Debug, Default, PartialEq, Clone, Copy)]
enum ScalingType { enum ScalingType {
Logarithmic, Logarithmic,
#[default] #[default]
Linear Linear,
} }
impl eframe::App for App { impl eframe::App for App {
@ -72,38 +75,65 @@ impl eframe::App for App {
egui::Grid::new("grid1") egui::Grid::new("grid1")
.min_col_width(125.0) .min_col_width(125.0)
.show(ui, |ui| { .show(ui, |ui| {
ui.label("calc interval start:").on_hover_text("the upper bound of the interval that functions are calculated in");; ui.label("calc interval start:").on_hover_text(
float_text_edit_singleline(ui, &mut self.a); "the lower bound of the interval that functions are calculated in",
ui.end_row(); );
float_text_edit_singleline(ui, &mut self.a);
ui.end_row();
ui.label("calc interval end:").on_hover_text("the upper bound of the interval that functions are calculated in"); ui.label("calc interval end:").on_hover_text(
float_text_edit_singleline(ui, &mut self.b); "the upper bound of the interval that functions are calculated in",
ui.end_row(); );
float_text_edit_singleline(ui, &mut self.b);
ui.label("N for omnes:").on_hover_text("the number of points used in calculation of omnes function"); ui.end_row();
int_text_edit_singleline(ui, &mut self.n_omnes);
ui.end_row();
ui.label("N for phi0:").on_hover_text("the number of points used in calculation of phi0 function"); ui.label("N for omnes:").on_hover_text(
int_text_edit_singleline(ui, &mut self.n_phi0); "the number of points used in calculation of omnes function",
ui.end_row(); );
int_text_edit_singleline(ui, &mut self.n_omnes);
ui.end_row();
ui.label("number of plotpoints:").on_hover_text("the number of points shown in the plot"); ui.label("N for phi0:").on_hover_text(
int_text_edit_singleline(ui, &mut self.n_points); "the number of points used in calculation of phi0 function",
ui.end_row(); );
}); int_text_edit_singleline(ui, &mut self.n_phi0);
ui.end_row();
ui.label("number of plotpoints:")
.on_hover_text("the number of points shown in the plot");
int_text_edit_singleline(ui, &mut self.n_points);
ui.end_row();
});
ui.separator(); ui.separator();
ui.vertical(|ui| { ui.vertical(|ui| {
ui.checkbox(&mut self.use_tan_omnes, "Use tan()-subst. for Omnes"); let omnes_cb =
ui.checkbox(&mut self.use_tan_phi0, "Use tan()-subst. for phi0"); ui.checkbox(&mut self.use_tan_omnes, "Use tan()-subst. for Omnes ('ot')");
ui.checkbox(&mut self.use_xsub_phi0, "Use second subst. for phi0"); if omnes_cb.clicked() {
egui::ComboBox::from_label("pick y-scaling").selected_text(format!("{:?}", self.scaling_type)).show_ui(ui, |ui| { self.calc_cache.omnes_lut = HashMap::with_hasher(Default::default()); // TODO split omnes_lut so that this is no longer necessary
ui.selectable_value(&mut self.scaling_type, ScalingType::Logarithmic, "log10"); }
ui.selectable_value(&mut self.scaling_type, ScalingType::Linear, "linear"); 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')");
ui.checkbox(
&mut self.use_reduced_integrand,
"Use reduced integrand for phi0",
);
egui::ComboBox::from_label("pick y-scaling")
.selected_text(format!("{:?}", self.scaling_type))
.show_ui(ui, |ui| {
ui.selectable_value(
&mut self.scaling_type,
ScalingType::Logarithmic,
"log10",
);
ui.selectable_value(
&mut self.scaling_type,
ScalingType::Linear,
"linear",
);
});
ui.horizontal(|ui| { ui.horizontal(|ui| {
let button_clear = ui.button("Clear canvas"); let button_clear = ui.button("Clear canvas");
if button_clear.clicked() { if button_clear.clicked() {
@ -111,99 +141,225 @@ impl eframe::App for App {
} }
egui::widgets::global_dark_light_mode_switch(ui); egui::widgets::global_dark_light_mode_switch(ui);
}); });
}); });
}); });
ui.separator(); ui.separator();
ui.label("You can calculate the following functions:"); ui.label("You can calculate the following functions:");
ui.horizontal( |ui| { ui.horizontal(|ui| {
for i in 0..self.plots_available.len() { for i in 0..self.plots_available.len() {
let button = ui.button( format!("{}", self.plots_available[i].clone())); let button = ui.button(format!("{}", self.plots_available[i].clone()));
if button.clicked() { if button.clicked() {
let x_values: Vec<f64> = (0..self.n_points).map(|x| -> f64 {f64::from(x) * ((self.b - self.a) / f64::from(self.n_points)) + self.a}).collect(); let x_values: Vec<f64> = (0..self.n_points)
.map(|x| -> f64 {
f64::from(x) * ((self.b - self.a) / f64::from(self.n_points))
+ self.a
})
.collect();
match i { match i {
0 => { 0 => {
let y_values_delta:Vec<f64> = if self.scaling_type == ScalingType::Linear { let y_values_delta: Vec<f64> = if self.scaling_type
x_values.iter().map(|&x| calc::delta_with_lut(&mut self.calc_cache, x)).collect() == ScalingType::Linear
{
x_values
.iter()
.map(|&x| calc::delta_with_lut(&mut self.calc_cache, x))
.collect()
} else { } else {
x_values.iter().map(|&x| calc::delta_with_lut(&mut self.calc_cache, x).log10()).collect() x_values
.iter()
.map(|&x| {
calc::delta_with_lut(&mut self.calc_cache, x)
.log10()
})
.collect()
}; };
self.plot_data.push(PlotCurve { self.plot_data.push(PlotCurve {
points: x_values.iter().zip(y_values_delta.iter()).map(|(&x,&y)| [x.powf(0.5),y]).collect(), points: x_values
name: "Delta".to_string() .iter()
.zip(y_values_delta.iter())
.map(|(&x, &y)| [x.powf(0.5), y])
.collect(),
name: format!("(#{}) Delta", self.plot_data.len() + 1),
}) })
} }
1 => { 1 => {
let y_values_omnes_integrand:Vec<f64> = if self.scaling_type == ScalingType::Linear { let y_values_omnes: Vec<f64> = if self.scaling_type
x_values.iter().map(|&x| calc::omnes_integrand_tan(&mut self.calc_cache, x, 100.0).abs()).collect() == ScalingType::Linear
{
print!("Calculating all plot points: ");
let bar =
indicatif::ProgressBar::new(u64::from(self.n_points));
let res = x_values
.iter()
.map(|&x| {
bar.inc(1);
calc::omnes(
&mut self.calc_cache,
x,
self.n_omnes,
self.use_tan_omnes,
)
.abs()
//.powi(2)
})
.collect();
bar.finish();
res
} else { } else {
x_values.iter().map(|&x| calc::omnes_integrand_tan(&mut self.calc_cache, x, 100.0).abs().log10()).collect() x_values
.iter()
.map(|&x| {
calc::omnes(
&mut self.calc_cache,
x,
self.n_omnes,
self.use_tan_omnes,
)
.abs()
.powi(2)
.log10()
})
.collect()
}; };
self.plot_data.push(PlotCurve { self.plot_data.push(PlotCurve {
points: x_values.iter().zip(y_values_omnes_integrand.iter()).map(|(&x,&y)| [x.powf(0.5),y]).collect(), points: x_values
name: "Omnes Integrand".to_string() .iter()
.zip(y_values_omnes.iter())
.map(|(&x, &y)| [x.powf(0.5), y])
.collect(),
name: format!(
"(#{}) Omnes, {} ot, n_o={}",
self.plot_data.len() + 1,
if self.use_tan_omnes { "w/" } else { "no " },
self.n_omnes
),
}) })
} }
2 => { 2 => {
let y_values_omnes:Vec<f64> = if self.scaling_type == ScalingType::Linear { let y_values_phi0: Vec<f64> =
x_values.iter().map(|&x| calc::omnes(&mut self.calc_cache, x, self.n_omnes, self.use_tan_omnes).abs().powi(2)).collect() if self.scaling_type == ScalingType::Linear {
} else { let bar =
x_values.iter().map(|&x| calc::omnes(&mut self.calc_cache, x, self.n_omnes, self.use_tan_omnes).abs().powi(2).log10()).collect() indicatif::ProgressBar::new(u64::from(self.n_points));
}; x_values
.iter()
.map(|&x| {
bar.inc(1);
calc::phi0(
&mut self.calc_cache,
x,
self.n_omnes,
self.n_phi0,
self.use_tan_omnes,
self.use_tan_phi0,
self.use_xsub_phi0,
self.use_reduced_integrand,
x == self.a,
)
})
.collect()
} else {
x_values
.iter()
.map(|&x| {
calc::phi0(
&mut self.calc_cache,
x,
self.n_omnes,
self.n_phi0,
self.use_tan_omnes,
self.use_tan_phi0,
self.use_xsub_phi0,
self.use_reduced_integrand,
x == self.a,
)
.log10()
})
.collect()
};
self.plot_data.push(PlotCurve { self.plot_data.push(PlotCurve {
points: x_values.iter().zip(y_values_omnes.iter()).map(|(&x,&y)| [x.powf(0.5),y]).collect(), points: x_values
name: "Omnes".to_string() .iter()
.zip(y_values_phi0.iter())
.map(|(&x, &y)| [x.powf(0.5), y])
.collect(),
name: format!(
"(#{}) phi0, {} ot, {} pt, {} ss, n_o={}, n_p={}",
self.plot_data.len() + 1,
if self.use_tan_omnes { "w/" } else { "no " },
if self.use_tan_phi0 { "w/" } else { "no " },
if self.use_xsub_phi0 { "w/" } else { "no " },
self.n_omnes,
self.n_phi0
),
}) })
} }
3 => { 3 => {
let y_values_phi0_integrand:Vec<f64> = if self.scaling_type == ScalingType::Linear { let y_values_phi0: Vec<f64> =
x_values.iter().map(|&x| calc::phi0_integrand_tan_xsub(&mut self.calc_cache, x, 10000.0, self.n_omnes, self.use_tan_omnes)).collect() if self.scaling_type == ScalingType::Linear {
} else { x_values
x_values.iter().map(|&x| calc::phi0_integrand_tan_xsub(&mut self.calc_cache, x, 10000.0, self.n_omnes, self.use_tan_omnes).log10()).collect() .iter()
}; .map(|&x| {
let val = calc::phi0(
&mut self.calc_cache,
x,
self.n_omnes,
self.n_phi0,
self.use_tan_omnes,
self.use_tan_phi0,
self.use_xsub_phi0,
self.use_reduced_integrand,
x == self.a,
) / calc::delta_with_lut(
&mut self.calc_cache,
x,
);
val
})
.collect()
} else {
x_values
.iter()
.map(|&x| {
let val = calc::phi0(
&mut self.calc_cache,
x,
self.n_omnes,
self.n_phi0,
self.use_tan_omnes,
self.use_tan_phi0,
self.use_xsub_phi0,
self.use_reduced_integrand,
x == self.a,
) / calc::delta_with_lut(
&mut self.calc_cache,
x,
);
val.log10()
})
.collect()
};
self.plot_data.push(PlotCurve { self.plot_data.push(PlotCurve {
points: x_values.iter().zip(y_values_phi0_integrand.iter()).map(|(&x,&y)| [x.powf(0.5),y]).collect(), points: x_values
name: "Phi0 Integrand".to_string() .iter()
.zip(y_values_phi0.iter())
.map(|(&x, &y)| [x.powf(0.5), y])
.collect(),
name: format!(
"(#{}) phi0 / delta, {} ot, {} pt, {} ss, n_o={}, n_p={}",
self.plot_data.len() + 1,
if self.use_tan_omnes { "w/" } else { "no " },
if self.use_tan_phi0 { "w/" } else { "no " },
if self.use_xsub_phi0 { "w/" } else { "no " },
self.n_omnes,
self.n_phi0
),
}) })
} }
4 => { 4_usize.. => {
let y_values_phi0: Vec<f64> = if self.scaling_type == ScalingType::Linear {
x_values.iter().map(|&x| calc::phi0(&mut self.calc_cache, x, self.n_omnes, self.n_phi0, self.use_tan_omnes, self.use_tan_phi0, self.use_xsub_phi0)).collect()
} else {
x_values.iter().map(|&x| calc::phi0(&mut self.calc_cache, x, self.n_omnes, self.n_phi0, self.use_tan_omnes, self.use_tan_phi0, self.use_xsub_phi0).log10()).collect()
};
self.plot_data.push(PlotCurve {
points: x_values.iter().zip(y_values_phi0.iter()).map(|(&x,&y)| [x.powf(0.5),y]).collect(),
name: "Phi0".to_string()
})
// let y_values_phi0 = calc::bit_pattern_randomness_tester(self.a, self.b, self.n_points as u64);
// self.plot_data.push(PlotCurve {
// points: (0..self.n_points).zip(y_values_phi0.iter()).map(|(x,&y)| [x as f64,y]).collect(),
// name: "Phi0".to_string()
// })
}
5 => {
let y_values_phi0: Vec<f64> = if self.scaling_type == ScalingType::Linear {
x_values.iter().map(|&x| {
let val = calc::phi0(&mut self.calc_cache, x.powi(2), self.n_omnes, self.n_phi0, self.use_tan_omnes, self.use_tan_phi0, self.use_xsub_phi0) / calc::delta_with_lut(&mut self.calc_cache, x.powi(2));
val
}).collect()
} else {
x_values.iter().map(|&x| {
let val = calc::phi0(&mut self.calc_cache, x.powi(2), self.n_omnes, self.n_phi0, self.use_tan_omnes, self.use_tan_phi0, self.use_xsub_phi0) / calc::delta_with_lut(&mut self.calc_cache, x.powi(2));
val.log10()
}).collect()
};
self.plot_data.push(PlotCurve {
points: x_values.iter().zip(y_values_phi0.iter()).map(|(&x,&y)| [x,y]).collect(),
name: "phi0 / delta".to_string()
})
}
6_usize.. => {
panic!("Not all buttons have been assigned data!"); panic!("Not all buttons have been assigned data!");
} }
} }
@ -211,10 +367,10 @@ impl eframe::App for App {
} }
}); });
}); });
if self.scaling_type != cur_scaling_type { if self.scaling_type != cur_scaling_type {
for curve in self.plot_data.iter_mut() { for curve in self.plot_data.iter_mut() {
for [_x,y] in curve.points.iter_mut() { for [_x, y] in curve.points.iter_mut() {
if self.scaling_type == ScalingType::Linear { if self.scaling_type == ScalingType::Linear {
*y = 10.0f64.powf(*y); *y = 10.0f64.powf(*y);
} else if self.scaling_type == ScalingType::Logarithmic { } else if self.scaling_type == ScalingType::Logarithmic {
@ -230,27 +386,42 @@ impl eframe::App for App {
.legend(Legend::default()); .legend(Legend::default());
let plot = match self.scaling_type { let plot = match self.scaling_type {
ScalingType::Linear => plot.custom_y_axes(vec![AxisHints::new_y().formatter(|mark, _l, _| scientific_notation(mark.value, false))]), ScalingType::Linear => plot.custom_y_axes(vec![AxisHints::new_y()
ScalingType::Logarithmic => {plot.y_grid_spacer(|input| { .formatter(|mark, _l, _| scientific_notation(mark.value, false))]),
let input_exp = GridInput {bounds: (inverse_log_map(input.bounds.0), inverse_log_map(input.bounds.1)), base_step_size: input.base_step_size}; ScalingType::Logarithmic => plot
.y_grid_spacer(|input| {
let input_exp = GridInput {
bounds: (
inverse_log_map(input.bounds.0),
inverse_log_map(input.bounds.1),
),
base_step_size: input.base_step_size,
};
let lin_gridmarks = (log_grid_spacer(10))(input_exp); let lin_gridmarks = (log_grid_spacer(10))(input_exp);
let mut log_gridmarks: Vec<GridMark> = Vec::new(); let mut log_gridmarks: Vec<GridMark> = Vec::new();
for v in lin_gridmarks.iter() { for v in lin_gridmarks.iter() {
if v.value - v.value.floor() < 0.9 { if v.value - v.value.floor() < 0.9 {
log_gridmarks.push(GridMark{value:log_map(v.value), step_size:v.step_size / 1.14}); log_gridmarks.push(GridMark {
} value: log_map(v.value),
step_size: v.step_size / 1.14,
});
}
} }
log_gridmarks log_gridmarks
}) })
.label_formatter(|_,point| format!("x = {:.5}\ny = {:.5}", point.x, 10.0f64.powf(point.y)).to_string()) .label_formatter(|_, point| {
.custom_y_axes(vec![AxisHints::new_y().formatter(|mark, _l, _| scientific_notation(mark.value, true))]) format!("x = {:.5}\ny = {:.5}", point.x, 10.0f64.powf(point.y)).to_string()
} })
}; .custom_y_axes(vec![AxisHints::new_y()
.formatter(|mark, _l, _| scientific_notation(mark.value, true))]),
};
plot.show(ui, |plot_ui| { plot.show(ui, |plot_ui| {
for i in 0..self.plot_data.len() { for i in 0..self.plot_data.len() {
plot_ui.line(Line::new(self.plot_data[i].points.clone()).name(self.plot_data[i].name.clone())); // is clone necessary? plot_ui.line(
Line::new(self.plot_data[i].points.clone())
.name(self.plot_data[i].name.clone()),
); // is clone necessary?
} }
}); });
}); });
@ -261,14 +432,23 @@ impl eframe::App for App {
fn main() { fn main() {
let mut app = App::default(); let mut app = App::default();
app.calc_cache = calc::Cache::from_file("./gauss_quad_lut.morello").unwrap(); app.calc_cache = calc::Cache::from_file("./gauss_quad_lut.morello").unwrap();
// let x_values: Vec<f64> = (0..app.n_points).map(|x| -> f64 {f64::from(x) * ((app.b - app.a) / f64::from(app.n_points)) + app.a}).collect(); // let x_values: Vec<f64> = (0..app.n_points)
// .map(|x| -> f64 { f64::from(x) * ((app.b - app.a) / f64::from(app.n_points)) + app.a })
// .collect();
// let t0 = Instant::now(); // let t0 = Instant::now();
// let y_values_phi0: Vec<f64> = x_values.iter().map(|&x| calc::phi0(&mut app.calc_cache, x, app.n_omnes, app.n_phi0)).collect(); // let y_values_phi0: Vec<f64> = x_values
// .iter()
// .map(|&x| calc::phi0(&mut app.calc_cache, x, app.n_omnes, app.n_phi0))
// .collect();
// std::hint::black_box(&y_values_phi0); // std::hint::black_box(&y_values_phi0);
// println!("{:?}", t0.elapsed()); // println!("{:?}", t0.elapsed());
eframe::run_native(
eframe::run_native("Morellus", eframe::NativeOptions::default(), Box::new(|_cc| Box::new(app))).unwrap(); "Titel",
eframe::NativeOptions::default(),
Box::new(|_cc| Box::new(app)),
)
.unwrap();
} }
#[cfg(target_arch = "wasm32")] #[cfg(target_arch = "wasm32")]
@ -277,10 +457,10 @@ fn main() {
wasm_bindgen_futures::spawn_local(async { wasm_bindgen_futures::spawn_local(async {
let res = gloo_net::http::RequestBuilder::new("/lut/gauss_quad_lut.morello") let res = gloo_net::http::RequestBuilder::new("/lut/gauss_quad_lut.morello")
.method(gloo_net::http::Method::POST) .method(gloo_net::http::Method::POST)
.send() .send()
.await .await
.unwrap(); .unwrap();
app.calc_cache = calc::Cache::from_slice(&res.binary().await.unwrap()).unwrap(); app.calc_cache = calc::Cache::from_slice(&res.binary().await.unwrap()).unwrap();
let web_options = eframe::WebOptions::default(); let web_options = eframe::WebOptions::default();
@ -296,22 +476,30 @@ fn main() {
}); });
} }
fn log_map(x:f64) -> f64 { fn log_map(x: f64) -> f64 {
x.floor() + (1.0 + (x - x.floor()) * 10.0).log10() x.floor() + (1.0 + (x - x.floor()) * 10.0).log10()
} }
fn inverse_log_map(x:f64) -> f64 { fn inverse_log_map(x: f64) -> f64 {
x.floor() + (10f64.powf(x - x.floor()) - 1.0) / 10.0 x.floor() + (10f64.powf(x - x.floor()) - 1.0) / 10.0
} }
fn scientific_notation(x:f64, log:bool) -> String { fn scientific_notation(x: f64, log: bool) -> String {
if log { if log {
format!("{}*10^{}", format!( format!(
"{:.8}", 10f64.powf(x - x.floor())) "{}*10^{}",
.trim_end_matches('0').trim_end_matches('.').to_string(), format!("{:.8}", 10f64.powf(x - x.floor()))
x.floor()).to_string() .trim_end_matches('0')
.trim_end_matches('.')
.to_string(),
x.floor()
)
.to_string()
} else { } else {
format!("{:.9}", x).trim_end_matches('0').trim_end_matches('.').to_string() format!("{:.9}", x)
.trim_end_matches('0')
.trim_end_matches('.')
.to_string()
} }
} }
@ -322,7 +510,6 @@ fn int_text_edit_singleline(ui: &mut egui::Ui, value: &mut u32) -> egui::Respons
*value = result; *value = result;
} else if tmp_value == "" { } else if tmp_value == "" {
*value = 0; *value = 0;
} }
res res
} }
@ -336,4 +523,4 @@ fn float_text_edit_singleline(ui: &mut egui::Ui, value: &mut f64) -> egui::Respo
*value = 0.0; *value = 0.0;
} }
res res
} }