make main single calculation

This commit is contained in:
Jan Bergen 2024-08-02 20:23:02 +02:00
parent d969de62f8
commit 1dbfafc999
3 changed files with 297 additions and 161 deletions

View File

@ -4,8 +4,6 @@ use std::collections::HashMap;
use std::io;
use std::rc::Rc;
use ahash::{AHasher, RandomState};
use std::{fs::File, io::BufReader};
use crate::gq_storage;
@ -217,7 +215,6 @@ pub fn omnes(cache: &mut Cache, s: f64, n: u32, use_tan: bool) -> Complex {
return *res;
}
}
if !cache.omnes_lut.contains_key(&n) {
cache
.omnes_lut
@ -275,6 +272,10 @@ pub fn phi0_integrand(
use_tan_phi0: bool,
use_xsub: bool,
_use_reduced_integrand: bool,
use_cutoff: bool,
cutoff_power: f64,
cutoff_s: f64,
cutoff_factor: Complex,
) -> f64 {
let sub = match (use_tan_phi0, use_xsub) {
(true, true) => s_tick.tan().powi(2) + S0,
@ -282,8 +283,18 @@ pub fn phi0_integrand(
(false, true) => s_tick.powi(2) + S0,
(false, false) => s_tick,
};
// let omnes_s = if use_cutoff && s >= cutoff_s {
// cutoff_factor * s.powf(cutoff_power)
// } else {
// omnes(cache, s, n_omnes, use_tan_omnes)
// };
let omnes_sub = if use_cutoff && sub >= cutoff_s {
cutoff_factor * sub.powf(cutoff_power)
// Complex::new((0.00002f64).powf(cutoff_power), 0.0)
} else {
omnes(cache, sub, n_omnes, use_tan_omnes)
};
let omnes_s = omnes(cache, s, n_omnes, use_tan_omnes);
let omnes_sub = omnes(cache, sub, n_omnes, use_tan_omnes);
let mut res = (omnes_sub / omnes_s).norm().ln() / (sub * (sub - s));
res *= match (use_tan_phi0, use_xsub) {
(true, true) => 2.0 / s_tick.cos().powi(2),
@ -304,6 +315,9 @@ pub fn phi0(
use_xsub: bool,
use_reduced_integrand: bool,
print_prog: bool,
use_cutoff: bool,
cutoff_power: f64,
cutoff_s: f64
) -> f64 {
if !cache.gauss_quad_lut.contains_key(&n_phi0) {
let gq_values = GaussLegendre::nodes_and_weights(n_phi0.try_into().unwrap());
@ -323,6 +337,9 @@ pub fn phi0(
.ln()
/ (s * (s - S0).powf(1.5))
};
let cutoff_factor = omnes(cache, cutoff_s, n_omnes, use_tan_omnes) * cutoff_s.powf(-cutoff_power);
let intgrl = integrate(
&*roots,
&*weights,
@ -344,6 +361,10 @@ pub fn phi0(
use_tan_phi0,
use_xsub,
use_reduced_integrand,
use_cutoff,
cutoff_power,
cutoff_s,
cutoff_factor,
) * s
/ PI
},
@ -352,8 +373,6 @@ pub fn phi0(
-(s - S0).sqrt() * (intgrl - analyt)
}
#[allow(unused)]
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];

View File

@ -7,7 +7,7 @@
use std::{
fs::File,
io::{self, BufReader, BufWriter, Read, Write, Seek, SeekFrom},
io::{self, BufReader, BufWriter, Read, Seek, SeekFrom, Write},
};
use gauss_quad::GaussLegendre;
@ -28,7 +28,10 @@ pub fn serialize(filepath: &str, params: Vec<(Vec<f64>, Vec<f64>)>) -> io::Resul
Ok(())
}
pub fn deserialize<R:Read + 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 values = Vec::new();
@ -41,7 +44,9 @@ pub fn deserialize<R:Read + Seek>(reader: &mut R, filter: impl Fn(u64) -> bool)
// skips stored floats for n that should not be stored
if !filter(n) {
reader.seek(SeekFrom::Current(i64::try_from(2*n*8).unwrap())).unwrap();
reader
.seek(SeekFrom::Current(i64::try_from(2 * n * 8).unwrap()))
.unwrap();
continue;
}

View File

@ -6,14 +6,14 @@ use indicatif;
use std::{collections::HashMap, time::Instant};
use egui_plot::{log_grid_spacer, AxisHints, GridInput, GridMark, Legend, Line, Plot};
use num::complex::ComplexFloat;
use num::{complex::ComplexFloat, Float};
mod calc;
mod gq_storage;
struct App {
plot_data: Vec<PlotCurve>,
plots_available: [String; 4],
plots_available: [String; 7],
calc_cache: calc::Cache,
scaling_type: ScalingType,
a: f64,
@ -26,8 +26,8 @@ struct App {
use_xsub_phi0: bool,
use_reduced_integrand: bool,
use_cutoff: bool,
cutoff_power: f64,
cutoff_s: f64,
cont_power: f64,
}
#[derive(Clone)]
@ -45,6 +45,9 @@ impl Default for App {
"Omnes".to_string(),
"Phi0".to_string(),
"Phi0 / delta".to_string(),
"s^cutoff_power".to_string(),
"phi0 integrand".to_string(),
"Phi0_cut / Phi0".to_string(),
],
calc_cache: calc::Cache::default(),
scaling_type: ScalingType::default(),
@ -62,8 +65,8 @@ impl Default for App {
// Cutoff params
use_cutoff: false,
cutoff_power: -1.0,
cutoff_s: 100.0,
cont_power: -1.0,
}
}
}
@ -82,42 +85,47 @@ impl eframe::App for App {
ui.vertical(|ui| {
ui.horizontal(|ui| {
egui::Grid::new("grid1")
.min_col_width(125.0)
.show(ui, |ui| {
ui.label("calc interval start:").on_hover_text(
"the lower bound of the interval that functions are calculated in",
);
float_text_edit_singleline(ui, &mut self.a, true);
ui.end_row();
ui.vertical(|ui| {
ui.label(egui::RichText::new("Boundary values"));
ui.add_space(10.0);
egui::Grid::new("grid1")
.min_col_width(125.0)
.show(ui, |ui| {
ui.label("calc interval start:").on_hover_text(
"the lower bound of the interval that functions are calculated in",
);
float_text_edit_singleline(ui, &mut self.a, true);
ui.end_row();
ui.label("calc interval end:").on_hover_text(
"the upper bound of the interval that functions are calculated in",
);
float_text_edit_singleline(ui, &mut self.b, true);
ui.end_row();
ui.label("calc interval end:").on_hover_text(
"the upper bound of the interval that functions are calculated in",
);
float_text_edit_singleline(ui, &mut self.b, true);
ui.end_row();
ui.label("N for omnes:").on_hover_text(
"the number of points used in calculation of omnes function",
);
int_text_edit_singleline(ui, &mut self.n_omnes, true);
ui.end_row();
ui.label("N for omnes:").on_hover_text(
"the number of points used in calculation of omnes function",
);
int_text_edit_singleline(ui, &mut self.n_omnes, true);
ui.end_row();
ui.label("N for phi0:").on_hover_text(
"the number of points used in calculation of phi0 function",
);
int_text_edit_singleline(ui, &mut self.n_phi0, true);
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, true);
ui.end_row();
});
ui.label("N for phi0:").on_hover_text(
"the number of points used in calculation of phi0 function",
);
int_text_edit_singleline(ui, &mut self.n_phi0, true);
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, true);
ui.end_row();
});
});
ui.separator();
ui.vertical(|ui| {
ui.label(egui::RichText::new("Use substitutions"));
ui.add_space(10.0);
let omnes_tan_cb =
ui.checkbox(&mut self.use_tan_omnes, "Use tan()-subst. for Omnes ('ot')");
if omnes_tan_cb.clicked() {
@ -155,18 +163,19 @@ impl eframe::App for App {
ui.separator();
ui.vertical(|ui| {
ui.label("Omnes Continuation with polynomial").on_hover_text(
ui.with_layout(egui::Layout::top_down(egui::Align::Min),|ui| {
ui.label(egui::RichText::new("Omnes Continuation with polynomial")).on_hover_text(
"If activated, from a user chosen cutoff value, the Omnes function will be replaced with s to a user chosen power, starting at the last Omnes value of the numerically calculated values",
);
ui.checkbox(&mut self.use_cutoff, "Use cutoff");
ui.checkbox(&mut self.use_cutoff, "Use cutoff for phi0");
ui.add_space(10.0);
egui::Grid::new("grid2").min_col_width(125.0).show(ui, |ui| {
ui.label("Power for continuation").on_hover_text("Omnes will be replaced with s^{this power}.");
float_text_edit_singleline(ui, &mut self.cont_power, self.use_cutoff);
ui.label("Power of cont. polynomial").on_hover_text("Omnes will be replaced with s^{this power}.");
float_text_edit_singleline(ui, &mut self.cutoff_power, true);
ui.end_row();
ui.label("Cutoff point").on_hover_text("Omnes will be replaced from this s-value onward.");
float_text_edit_singleline(ui, &mut self.cutoff_s, self.use_cutoff);
float_text_edit_singleline(ui, &mut self.cutoff_s, true);
ui.end_row();
});
})
@ -268,49 +277,41 @@ impl eframe::App for App {
})
}
2 => {
let y_values_phi0: Vec<f64> =
if self.scaling_type == ScalingType::Linear {
let bar =
indicatif::ProgressBar::new(u64::from(self.n_points));
let t0 = Instant::now();
let values = 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();
println!("{:?}", t0.elapsed());
values
} 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()
};
let y_values_phi0: Vec<f64> = {
let bar =
indicatif::ProgressBar::new(u64::from(self.n_points));
let t0 = Instant::now();
let values = x_values
.iter()
.map(|&x| {
bar.inc(1);
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,
self.use_cutoff,
self.cutoff_power,
self.cutoff_s,
);
match self.scaling_type {
ScalingType::Linear => {
val
}
ScalingType::Logarithmic => {
val.log10()
}
}
}
).collect();
println!("{:?}", t0.elapsed());
values
};
self.plot_data.push(PlotCurve {
points: x_values
.iter()
@ -326,52 +327,40 @@ impl eframe::App for App {
self.n_omnes,
self.n_phi0
),
})
});
}
3 => {
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,
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()
let y_values_phi0: Vec<f64> = {
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,
self.use_cutoff,
self.cutoff_power,
self.cutoff_s,
) / calc::delta_with_lut(
&mut self.calc_cache,
x,
);
match self.scaling_type {
ScalingType::Linear => {
val
}
ScalingType::Logarithmic => {
val.log10()
}
}
})
.collect()
};
self.plot_data.push(PlotCurve {
points: x_values
@ -388,10 +377,130 @@ impl eframe::App for App {
self.n_omnes,
self.n_phi0
),
});
}
4 => {
let cutoff_factor = calc::omnes(&mut self.calc_cache, self.cutoff_s, self.n_omnes, self.use_tan_omnes).abs() * self.cutoff_s;
println!("{:?}", self.scaling_type);
let y_values_func: Vec<f64> = x_values.iter().map(|&x| {
let val = (cutoff_factor*x.powf(self.cutoff_power)).powi(2);
match self.scaling_type {
ScalingType::Linear => {
val
}
ScalingType::Logarithmic => {
val.log10()
}
}
}).collect();
self.plot_data.push(PlotCurve {
points: x_values
.iter()
.zip(y_values_func.iter())
.map(|(&x, &y)| [x.powf(0.5), y])
.collect(),
name:"cont. polynomial".to_string()
});
}
5 => {
let cutoff_factor = calc::omnes(&mut self.calc_cache, self.cutoff_s, self.n_omnes, self.use_tan_omnes) * self.cutoff_s;
println!("{:?}", cutoff_factor);
let y_values_func: Vec<f64> = x_values
.iter()
.map(|&x| {
let val = calc::phi0_integrand(
&mut self.calc_cache,
x,
2.0,
self.n_omnes,
self.use_tan_omnes,
self.use_tan_phi0,
self.use_xsub_phi0,
self.use_reduced_integrand,
self.use_cutoff,
self.cutoff_power,
self.cutoff_s,
cutoff_factor);
match self.scaling_type {
ScalingType::Linear => {
val
}
ScalingType::Logarithmic => {
val.log10()
}
}
}
).collect();
self.plot_data.push(PlotCurve {
points: x_values
.iter()
.zip(y_values_func.iter())
.map(|(&x, &y)| [x.powf(0.5), y])
.collect(),
name: format!(
"(#{}) phi0 integrand",
self.plot_data.len() + 1,
),
})
}
4_usize.. => {
panic!("Not all buttons have been assigned data!");
6 => {
let bar = indicatif::ProgressBar::new(u64::from(self.n_points));
let y_values: Vec<f64> =
x_values.iter().map(|&x| {
bar.inc(1);
let phi0_nocut = calc::phi0(
&mut self.calc_cache,
3.0,
self.n_omnes,
self.n_phi0,
self.use_tan_omnes,
self.use_tan_phi0,
self.use_xsub_phi0,
self.use_reduced_integrand,
false,
false,
self.cutoff_power,
x
);
let phi0_cut = calc::phi0(
&mut self.calc_cache,
3.0,
self.n_omnes,
self.n_phi0,
self.use_tan_omnes,
self.use_tan_phi0,
self.use_xsub_phi0,
self.use_reduced_integrand,
false,
true,
self.cutoff_power,
x
);
match self.scaling_type {
ScalingType::Linear => {
(phi0_cut / phi0_nocut - 1.0).abs()
}
ScalingType::Logarithmic => {
(phi0_cut / phi0_nocut - 1.0).abs().log10()
}
}
}).collect();
self.plot_data.push(PlotCurve {
points: x_values
.iter()
.zip(y_values.iter())
.map(|(&x, &y)| [x.powf(0.5), y])
.collect(),
name: format!(
"(#{}) phi0_cut / phi0, cutoff power: {}",
self.plot_data.len() + 1,
self.cutoff_power,
),
});
}
7_usize.. => {
panic!("This button has no assigned purpose!");
}
}
}
@ -464,30 +573,33 @@ fn main() {
let mut app = App::default();
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 bar = indicatif::ProgressBar::new(u64::from(app.n_points));
let values: Vec<f64> = x_values
.iter()
.map(|&x| {
bar.inc(1);
calc::phi0(
&mut app.calc_cache,
x,
1000,
1000,
app.use_tan_omnes,
app.use_tan_phi0,
app.use_xsub_phi0,
app.use_reduced_integrand,
x == 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 bar = indicatif::ProgressBar::new(u64::from(app.n_points));
// let values: Vec<f64> = x_values
// .iter()
// .map(|&x| {
// bar.inc(1);
// calc::phi0(
// &mut app.calc_cache,
// x,
// 1000,
// 1000,
// app.use_tan_omnes,
// app.use_tan_phi0,
// app.use_xsub_phi0,
// app.use_reduced_integrand,
// x == app.a,
// )
// })
// .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);
// eframe::run_native(
// "Titel",
// "Omnes Calculator",
// eframe::NativeOptions::default(),
// Box::new(|_cc| Box::new(app)),
// )